C C multiparticle-scattering Fortran code gmm02s.f C calculates radiative scattering by an arbitrary external aggregate of C homogeneous and/or core-mantle spheres C the usage is exactly the same as codes gmm01f.f, gmm01s.f, and C gmm02f.f C this code is essentially the same as gmm02f.f except that it is less C computer-memory demanding and slower than gmm02f.f C Yu-lin Xu C released to public 9/2001 C available online at http://www.astro.ufl.edu/~xu C last updated 5/2003 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----------------------------------------------------------------------- PROGRAM gmm02s 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 component sphere in an aggregate C an example of "gmm01f.par": parameter (nLp=10000,np=20) C---------------------------------------------------------------------- parameter (nmp=np*(np+2),nmp0=(np+1)*(np+4)/2) parameter (nangmax=181,MOR=181,ncmax=360) parameter (ni0=np*(np+1)*(2*np+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) integer u,v,u0,nmax(nLp),uvmax(nLp),ind(nLp) double precision k,r0(9,nLp),x(nLp),xc(nLp), + dang(nangmax),r00(3,nLp),betar(MOR), + thetr(MOR),phair(MOR),smue(4,4),mue(4,4,ncmax,nangmax), + besj(0:2*np+1),besy(0:2*np+1),i11(nangmax), + i21(nangmax),i22(nangmax),i12(nangmax),inat(nangmax), + pol(nangmax),cscaxi(nLp),cscayi(nLp),cextxi(nLp), + cextyi(nLp),cabsxi(nLp),cabsyi(nLp),cexti(nLp),cabsi(nLp), + cscai(nLp),assymi(nLp),assymxi(nLp),assymyi(nLp), + cprxi(nLp),cpryi(nLp),cpri(nLp),drot(nrc), + c0i(nLp),c1i(nLp) complex*16 A,B,cmz,Aj,Bj,A2,B2,Aj2,Bj2,A0,B0,ci,cin, + atr(2,np,nmp),at(nmp),bt(nmp), + ek(np),p0(nLp,nmp),ref(nLp),refc(nLp), + q0(nLp,nmp),an(np),bn(np),aMie(nLp,np),bMie(nLp,np), + as(nLp,nmp),bs(nLp,nmp),as0(nLp,nmp),bs0(nLp,nmp), + asc(nLp,nmp),bsc(nLp,nmp),as1(nLp,nmp),bs1(nLp,nmp), + ast(nLp,nmp),bst(nLp,nmp),asp(nLp,nmp),bsp(nLp,nmp), + asv(nLp,nmp),bsv(nLp,nmp), + s2x(ncmax,nangmax),s4x(ncmax,nangmax), + s3y(ncmax,nangmax),s1y(ncmax,nangmax),atj(nmp),btj(nmp) CHARACTER FLNAME*20,fileout*20,fileout1*19,fileout2*21, + 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:2*(np+2)) 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) gcs=0.d0 gcv=0.d0 idpq=0 OPEN(UNIT=1,FILE='gmm01f.in',STATUS='OLD') READ(1,'(a20)') FLNAME C----------------------------------------------------------------------- C FLNAME - the input file name for the aggregate configuration C the first line in this file is the incident wavelength, C the second line is the number of spheres, C rest lines provide coordinates of sphere-centers, radii, C and refractive indexes of mantle and core; each line C contains nine real numbers: C x, y, z, r, Re(m), Im(m), rc, Rec(m), Imc(m) C rc=0.0 for a homogeneous sphere C----------------------------------------------------------------------- write(6,'(a12,a20)') 'input file: ',FLNAME READ(1,*) nbeta,nthet,nphai C----------------------------------------------------------------------- C The product of nbeta, nthet, and nphai is the total number of C orientations to be averaged for the input sphere-aggregate. C In the Cartesian coordinate system, the direction of propagation of C the incident plane wave defines the positive z-axis, the x-z plane is C the scattering plane. C nbeta - # of rotation around the z-axis, C this rotation angle corresponds to the Euler angle of "alpha" C nthet - # of rotation around the y-axis, C corresponding to the Euler angle of "beta" C nphai - # of rotation by the Euler angle of "gamma" C The definitions of the three Euler angles follow the convention used C by Edmonds ["Angular Momentum in Quantum Mechanics," pp.6-8 (1957)]. C examples: C When only a single orientation needs to be calculated: C nbeta=nthet=nphai=1 (the input line will be 1,1,1 or 1 1 1) C When an aggregate rotates around the y-axis (i.e., rotates in the C scattering plane), the number of orientations of the aggregate to be C calculated is 19, then nbeta=1, nthet=19, nphai=1 (1,19,1 or 1 19 1) C C *** In addition to the output file of mueller.out for the scattering C Mueller matrix, there is also an output file of amp.out for the C amplitude scattering matrix when nbeta=nthet=nphai=1 (i.e., when C a single fixed-orientation is calculated). C----------------------------------------------------------------------- write(6,'(a19,3i5)') 'nbeta,nthet,nphai: ', + nbeta,nthet,nphai if(nbeta.gt.MOR.or.nthet.gt.MOR.or.nphai.gt.MOR) then write(6,*) '*** parameter MOR too small ***' write(6,*) 'MOR must be >nbeta,nthet,nphai given above' write(6,*) 'Please change MOR in the parameter line of' write(6,*) 'the main code, recompile, then try again' stop endif if(nbeta*nphai.gt.1) then idran=1 else idran=0 endif nram=nbeta*nthet*nphai if(nram.lt.1) then write(6,*) + 'please check (nbeta,nthet,nphai) in gmm01f.in' stop endif READ(1,*) betami,betamx,thetmi,thetmx,phaimi,phaimx C----------------------------------------------------------------------- C betami,betamx,thetmi,thetmx,phaimi,phaimx are in degrees C betami,betamx -- the range of the Euler angle of rotation (alpha) C to be calculated C thetmi,thetmx -- the range of the Euler angle of rotation (beta) C to be calculated C phaimi,phaimx -- the range of the Euler angle of rotation (gamma) C to be calculated C an example: C When only a fixed orientation needs to be calculated, betami=betamx, C thetmi=thetmx, and phaimi=phaimx, this line could be, e.g., C (0. 0. 0. 0. 0. 0.), (0. 0. 90. 90. 0. 0.), (30. 30. 40. 40. 45. 45.) C------------------------------------------------------------------------ if(nbeta.eq.1) betamx=betami if(nthet.eq.1) thetmx=thetmi if(nphai.eq.1) phaimx=phaimi write(6,'(a24,6f7.2)') 'Ranges of Euler angles: ', + betami,betamx,thetmi,thetmx,phaimi,phaimx READ(1,*) idMie C------------------------------------------------------------------------ C idMie=1: calculating only coherent Mie-scattering, no interaction C------------------------------------------------------------------------ write(6,'(a7,i3)') 'idMie: ',idMie READ(1,*) idd C------------------------------------------------------------------------ C When calculating a single orientation, put idd=0 C When idd=1, the number of orientations to be calculated is doubled. C Each orientation is coupled with the orientation that the aggregate C is rotated by 90 degrees around the z axis. This insures that the C averaged polarizations are zero when the scattering angle is either 0 C or 180 degrees, as is suppossed to be for an average over random C orientations. C------------------------------------------------------------------------ if(nram.eq.1) idd=0 READ(1,*) idc,iseed C------------------------------------------------------------------------ C idc for choosing in the three schemes for an orientation average C idc=1 to devide the range of [thetmi,thetmx] by the cos(theta) scheme C (theta is the scattering angle, corresponding to the Euler angle C of ritation "beta") C idc=0 to devide the range of [thetmi,thetmx] by "degrees" C idc=-1 to pick up all the three Euler angles of rotation randomly using C a random number generator C iseed is the seed number for the random number generator, used only C when idc=-1 (i.e., when idc=1 or 0, iseed has no function) C iseed can be an arbitrary positive integer C when calculating only one single orientation, i.e., C nram=nbeta*nthet*nphai=1, idc and iseed have no function C------------------------------------------------------------------------ if(idpq.eq.1) then nram=nphai+nthet idc=0 idd=0 nbeta=1 betami=0.d0 betamx=betami thetr(nthet+1)=0.d0 phair(nphai+1)=0.d0 endif write(6,'(a5,i3)') 'idd: ',idd if(idd.eq.1) nram=2*nram write(6,'(a36,i5)') '# of orientations to be averaged: ',nram if(idc.lt.0) then write(6,'(a11,i4,i12)') 'idc,iseed: ',idc,iseed else write(6,'(a5,i3)') 'idc: ',idc endif READ(1,*) factor1,factor2,MXINT C------------------------------------------------------------------------ C factor1 and factor2 are numerical factors used for improving the C convergence behavior of the iterative solution process in solving C interacting equations, which are in the range of [0,1] C factor1 is for x-polarized incident plane wave and factor2 for C y-polarized incident wave C MXINT is the maximum number of iterations allowed in the iterative C solution process C This code uses two alternative methods in the iterative solution of C interacting equations: iteration scheme [see Fuller and Kattawar, Opt. C Lett. 13, 90 (1988); Xu, Appl. Opt. 34, 4573 (1995)] and BI-CGSTAB, the C stabilized Bi-Conjugate Gradient [see H.A. van der Vorst, SIAM J. Sci. C Stat. Comput. 13, 631, (1992)]. C When factor1=0 or factor2=0, the code directly goes to BI-CGSTAB C without using the other iteration scheme. C When factor1=factor2=1 it is equivalent to Fuller and Kattawar's C order-of-scattering method. C In many cases, a divergence will occur when factor1=factor2=1. Then, C setting factor1,factor2<1 may help to converge to a numerical solution. C When MXINT is exceeded, it will automatically switch to BI-CGSTAB. C C ******** In most circumstances the BI-CGSTAB is much more efficient C than the other iterative schemes so that factor1=factor2=0 is the best C choice for the iterative solution of the boundary condition equations. C When a user would like to try other iterative schemes, just simply C comment out the following two lines factor1=0 and factor2=0 and C the user-specified input values of factor1 and factor2 will be used . C 10/20/02 C------------------------------------------------------------------------ factor1=0 factor2=0 write(6,'(a,2f6.2)') 'Numerical factors for convergence:', + factor1,factor2 write(6,'(a37,i5)') 'Maximum iterations allowed:',MXINT READ(1,*) NADD C------------------------------------------------------------------------ C NADD is the number of terms to add to the scattering orders required C by the Wiscombe's criterion, which can be negative or positive in the C range of [-9,99] C Normally, set NADD=0 C------------------------------------------------------------------------ write(6,'(a41,i3)') + 'Scat. orders added to Wiscombe criterion:',NADD READ(1,*) eps,small if(eps.gt.1.d-20) eps=1.d-20 C------------------------------------------------------------------------ C eps: error tolerance for determining single-sphere field-expansion C truncation, default: 1.d-20 C (the default value of 1.d-20 allows to use Wiscombi's criterion) C small: error tolerance for the iterative solution process for solving C the interacting scattering coefficients (1.d-6) C------------------------------------------------------------------------ write(6,'(a35,e10.2)') + 'error tolerance for Mie-expansions:',eps write(6,'(a22,e10.2)') 'Convergence criterion:',small READ(1,*) fint C------------------------------------------------------------------------ C fint is the interaction index in the range of [0,1] (default: 0.02) 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 f0, the number of azimuth angles in the range of C (0,360) degrees that will be calculated in addition to 0 (360) C degrees is npng-1 with npng=360/pang. For example, when C pang=180, npng=2, the azimuth angles 0 (360) and 180 degrees C will be calculated. In the output for the Mueller matrix at C each azimuth angle, there are nang2 (=2*nang-1) sets of the C 16 elements. C------------------------------------------------------------------------ write(6,'(a32,f7.3)') + 'scat.-angle-interval in output: ',sang if(sang.le.0.d0) sang=1.d0 nang=90.d0/sang+1.d0 nang2=2*nang-1 if(nang2.gt.nangmax) then write(6,*) 'sang too small' write(6,*) + 'please increase sang in the input file gmm01f.in' write(6,*) ' and try again, or' write(6,*) + ' increase nangmax in the parameter line of the' write(6,*) ' main code, recompile, then try again' stop endif write(6,'(a,f8.3)') + 'azimuth-angle-interval in Mueller matrix output: ',pang if(pang.lt.0.0001d0) then npng=1 else npng=360.0d0/pang endif if(npng.gt.ncmax) then write(6,*) 'pang too small' write(6,*) + 'please increase pang in the input file gmm01f.in' write(6,*) + 'and try again, or increase ncmax in the parameter' write(6,*) + ' line of the main code, recompile, then try again' stop endif close(1) write(6,'(/)') OPEN(UNIT=2,FILE=FLNAME,STATUS='OLD') READ(2,*) w C------------------------------------------------------------------------ C w -- incident wavelength C------------------------------------------------------------------------ READ(2,*) nL C------------------------------------------------------------------------ C nL -- number of spheres in the aggregate C------------------------------------------------------------------------ if(nL.gt.nLp) then write(6,*) 'Parameter nLp too small, must be >', nL write(6,*) + 'Change nLp in gmm01f.par, recompile, then try again' stop endif if(nL.eq.1) then idMie=1 betami=0.d0 betamx=0.d0 thetmi=0.d0 thetmx=0.d0 phaimi=0.d0 phaimx=0.d0 nbeta=1 nthet=1 nphai=1 nram=1 idc=0 idd=0 endif C------------------------------------------------------------------------ C input the configuration and particle parameters for the aggregate C each line includes 9 numbers: C r0(1,i),r0(2,i),r0(3,i): C x-, y-, z-coordinates of the sphere-center i C r0(4,i),r0(5,i),r0(6,i): C the radius of the mantle of sphere i in the same unit of the C incident wavelength, the real and imaginary parts of the refractive C index of the mantle C r0(7,i),r0(8,i),r0(9,i): C the radius of the core of sphere i in the same unit of the C incident wavelength, the real and imaginary parts of the refractive C index of the core C------------------------------------------------------------------------ do 1 i=1,nL read(2,*,err=10) (r0(j,i),j=1,9) x0=r0(1,i) y0=r0(2,i) z0=r0(3,i) r00(1,i)=x0 r00(2,i)=y0 r00(3,i)=z0 r0(6,i)=dabs(r0(6,i)) r0(9,i)=dabs(r0(9,i)) ref(i)=dcmplx(r0(5,i),r0(6,i)) refc(i)=dcmplx(r0(8,i),r0(9,i)) 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,'(/)') betami=betami*pih/90.d0 betamx=betamx*pih/90.d0 thetmi=thetmi*pih/90.d0 thetmx=thetmx*pih/90.d0 phaimi=phaimi*pih/90.d0 phaimx=phaimx*pih/90.d0 if(idc.gt.0) then call orientcd(betami,betamx,thetmi,thetmx,phaimi,phaimx, + MOR,MOR,MOR,nbeta,nthet,nphai,betar,thetr,phair) else call orientud(betami,betamx,thetmi,thetmx,phaimi,phaimx, + MOR,MOR,MOR,nbeta,nthet,nphai,betar,thetr,phair) endif fileout='gmm02s.out' fileout1=fileout fileout2='pq'//fileout1 if(idMie.eq.1) then write(6,*) + '*** Calculating only coherent Mie-scattering ***' write(6,*) + '*** No interaction included ********************' endif do i=1,nL x(i)=k*r0(4,i) xc(i)=k*r0(7,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 if(xc(i).eq.xc(i-1).and.refc(i).eq.refc(i-1)) then nmax(i)=nmax(i-1) uvmax(i)=uvmax(i-1) goto 15 endif endif 12 write(6,'(a,i3,a,f10.4)') 'sphere #',i, + ' individual size parameter: ',x(i) c c calculating scattering coefficients for each component spherical c particles c subroutine scoatabd.f used here is originally SCSERD.FOR written by c R.T. Wang and W.X. Wang [see R.T. Wang and W.X. Wang, "Scattering by C Arbitrarily Large Homogeneous/Concentric Speres - Exact Theory with C Use of New Efficient Algorithms," in Proc. of the 1985 CRDC C Scientific Conference on Obscuration and Aerosol Research, R. Kohl, C ed. (Army CRDEC-SP-86019, Aberdeen, MD 1986), pp. 381-409], which C uses ratio method of Wang and van der Hulst in calculating C Riccati-Bessel functions in Mie-calculations [see Wang and van der C Hulst, Appl. Opt. 30, 106 (1991), Xu, Gustafson, Giovane, Blum, and c Tehranian, Phys. Rev. E 60, 2347 (1999)] c ratio=xc(i)/x(i) if(ratio.gt.1.d0) then write(6,*) 'size of core >mantle for particle ',i stop endif if(ratio.lt.1.d-14) then ratio=0.d0 endif call scoatabd(x(i),ratio,r0(8,i),-r0(9,i),r0(5,i), + -r0(6,i),np,an,bn,NADD,nmax(i)) 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 cextx=0.d0 cexty=0.d0 cabsx=0.d0 cabsy=0.d0 cscax=0.d0 cscay=0.d0 cprx=0.d0 cpry=0.d0 cbakx=0.d0 cbaky=0.d0 do i=1,nL cextxi(i)=0.d0 cextyi(i)=0.d0 cabsxi(i)=0.d0 cabsyi(i)=0.d0 cscaxi(i)=0.d0 cscayi(i)=0.d0 cprxi(i)=0.d0 cpryi(i)=0.d0 enddo do i=1,nang2 i11(i)=0.d0 i21(i)=0.d0 i22(i)=0.d0 i12(i)=0.d0 do jc=1,npng do j=1,4 do m=1,4 mue(j,m,jc,i)=0.d0 enddo enddo enddo enddo iram=0 if(idpq.eq.1) then open(12,file=fileout2,status='unknown') write(12,*) 'input file: ',FLNAME endif write(6,'(/)') write(6,*) 'original 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) if(idpq.eq.1) then nphaic=nphai+1 phair(nphaic)=0.d0 else nphaic=nphai endif do ibeta=1,nbeta do iphai=1,nphaic if(idpq.eq.1.and.iphai.lt.nphaic) then nthetc=1 else nthetc=nthet endif do ithet=1,nthetc if(idc.lt.0) then betar(ibeta)=(betamx-betami)*ran1d(iseed) phair(iphai)=(phaimx-phaimi)*ran1d(iseed) thetr(ithet)=(thetmx-thetmi)*ran1d(iseed) endif do 19 irot=1,2 if(idd.ne.1.and.irot.eq.2) goto 19 iram=iram+1 if(irot.eq.1) then alph=0.d0 else alph=pih endif ca=dcos(alph) sa=dsin(alph) beta=betar(ibeta) cb=dcos(beta) sb=dsin(beta) do i=1,nL x0=r00(1,i) y0=r00(2,i) r0(1,i)=cb*x0-sb*y0 r0(2,i)=sb*x0+cb*y0 enddo phai=phair(iphai) thet=thetr(ithet) if(idpq.eq.1.and.nthetc.eq.1) thet=0.d0 cb=dcos(phai) sb=dsin(phai) cz=dcos(thet) sz=dsin(thet) if(iram.eq.1.or.iram/50*50.eq.iram) then write(6,'(a,2i5)') 'iram & nram: ',iram,nram endif do i=1,nL x0=r0(1,i) y0=r0(2,i) z0=r00(3,i) r0(1,i)=ca*cz*x0-(ca*sz*sb+sa*cb)*y0 + +(ca*sz*cb-sa*sb)*z0 r0(2,i)=sa*cz*x0-(sa*sz*sb-ca*cb)*y0 + +(sa*sz*cb+ca*sb)*z0 r0(3,i)=-sz*x0-cz*sb*y0+cz*cb*z0 enddo if(iram.eq.1.or.iram/50*50.eq.iram) then i=1 write(6,'(i5,3f14.5)') + i,r0(1,i),r0(2,i),r0(3,i) i=nL write(6,'(i5,3f14.5)') + i,r0(1,i),r0(2,i),r0(3,i) endif C------------------------------------------------------------------------ C calculating constants, Gaunt, rotational and translation coefficients C------------------------------------------------------------------------ n0=nmax0+2 fnr(0)=0.d0 do n=1,2*n0 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) indpol=0 factor=factor1 if(iram/50*50.eq.iram) then write(6,'(a8,i4,a32)') + ' orien.#',iram,' Solving for x-pol. inci. state' endif 18 do imn=1,nmp do i=1,nL p0(i,imn)=0.d0 q0(i,imn)=0.d0 as(i,imn)=0.d0 bs(i,imn)=0.d0 enddo enddo C-------------------------------------------------- C calculating incident wave expansion coefficients C-------------------------------------------------- do 20 i=1,nL cz=dcos(k*r0(3,i)) sz=dsin(k*r0(3,i)) cmz=0.5d0*dcmplx(cz,sz) do 21 n=1,nmax(i) imn=n*n+n+1 A=fnr(2*n+1)*cmz p0(i,imn)=aMie(i,n)*A q0(i,imn)=bMie(i,n)*A p0(i,imn-2)=-p0(i,imn) q0(i,imn-2)=q0(i,imn) if(indpol.gt.1) then p0(i,imn)=p0(i,imn)*cin q0(i,imn)=q0(i,imn)*cin p0(i,imn-2)=p0(i,imn) q0(i,imn-2)=-q0(i,imn) endif as(i,imn)=p0(i,imn) bs(i,imn)=q0(i,imn) as(i,imn-2)=p0(i,imn-2) bs(i,imn-2)=q0(i,imn-2) 21 continue 20 continue C-------------------------------------------------------------- c begins iteration process to solve the interaction equations c for partial interacting scatteing coefficients c factor1,factor2=0: BI-CGSTAB [see H.A. van der Vorst, SIAM, c J. Sci. Stat. Comput. 13, 631 (1992)] c factor1=factor2=1: order-of-scattering [Fuller and Kattawar, c Opt. Lett. 13, 90 & 1063 (1988)] c 0' write(6,'(2x,6e13.5)') cext,cabs,csca,cbak,cext-cpr,assym cscax=0.d0 cscay=0.d0 cextx=0.d0 cexty=0.d0 c cabsx=0.d0 c cabsy=0.d0 cprx=0.d0 cpry=0.d0 do i=1,nL cscax=cscax+cscaxi(i) cscay=cscay+cscayi(i) cextx=cextx+cextxi(i) c cabsx=cabsx+cabsxi(i) cexty=cexty+cextyi(i) c cabsy=cabsy+cabsyi(i) cprx=cprx+cprxi(i) cpry=cpry+cpryi(i) enddo assymx=cprx/cscax assymy=cpry/cscay assym0=0.5d0*(assymx+assymy) cscax=2.0d0*twopi*cscax/cz cscay=2.0d0*twopi*cscay/cz csca=0.5d0*(cscax+cscay) cextx=twopi*cextx/cz cexty=twopi*cexty/cz cext=0.5d0*(cextx+cexty) cabsx=2.0d0*twopi*cabsx/cz cabsy=2.0d0*twopi*cabsy/cz cabs=0.5d0*(cabsx+cabsy) cabsx=cextx-cscax cabsy=cexty-cscay cabs=cext-csca cprx=2.0d0*twopi*cprx/cz cpry=2.0d0*twopi*cpry/cz cpr=0.5d0*(cprx+cpry) assym=cpr/csca write(6,'(2x,6e13.5)') cext,cabs,csca,cbak,cext-cpr,assym do i=1,nL c cabsxi(i)=4.0d0*pione*cabsxi(i)/cz c cabsyi(i)=4.0d0*pione*cabsyi(i)/cz cextxi(i)=2.0d0*pione*cextxi(i)/cz cextyi(i)=2.0d0*pione*cextyi(i)/cz cscaxi(i)=4.0d0*pione*cscaxi(i)/cz cscayi(i)=4.0d0*pione*cscayi(i)/cz cabsxi(i)=cextxi(i)-cscaxi(i) cabsyi(i)=cextyi(i)-cscayi(i) cprxi(i)=4.0d0*pione*cprxi(i)/cz cpryi(i)=4.0d0*pione*cpryi(i)/cz cscai(i)=0.5d0*(cscaxi(i)+cscayi(i)) cexti(i)=0.5d0*(cextxi(i)+cextyi(i)) cabsi(i)=0.5d0*(cabsxi(i)+cabsyi(i)) cpri(i)=0.5d0*(cprxi(i)+cpryi(i)) cpri(i)=cscai(i)+cabsi(i)-cpri(i) assymxi(i)=cprxi(i)/cscaxi(i) assymyi(i)=cpryi(i)/cscayi(i) assymi(i)=0.5d0*(cprxi(i)+cpryi(i))/csca cprxi(i)=cscaxi(i)+cabsxi(i)-cprxi(i) cpryi(i)=cscayi(i)+cabsyi(i)-cpryi(i) write(6,'(i5,5e15.6)') i,cexti(i),cabsi(i), + cscai(i),cpri(i),assymi(i) enddo write(6,'(/,a)') 'efficiencies for radiation pressure' write(6,'(5x,6e13.5)') assym,assym0 write(6,'(5x,6e13.5)') assym,cext-cpr,assymx,cextx-cprx, + assymy,cexty-cpry do i=1,nL write(6,'(i5,6e13.5)') i,assymi(i),cpri(i), + assymxi(i),cprxi(i),assymyi(i),cpryi(i) enddo betami=betami*90.d0/pih betamx=betamx*90.d0/pih thetmi=thetmi*90.d0/pih thetmx=thetmx*90.d0/pih phaimi=phaimi*90.d0/pih phaimx=phaimx*90.d0/pih flout='cr'//fileout open(33,file=flout,status='unknown') write(33,'(a20,a47)') + flout,'(Total and individual-particle cross sections)' write(33,'(a32,2x,a22)') + 'input sphere-aggregate filename:',FLNAME write(33,'(a19,3i5)') 'nbeta,nthet,nphai: ', + nbeta,nthet,nphai write(33,'(a24,6f7.2)') 'Ranges of Euler angles: ', + betami,betamx,thetmi,thetmx,phaimi,phaimx write(33,'(a28,i5)') '# of orientations averaged: ',nram write(33,'(12x,a4,11x,a4,11x,a4,11x,a3,8x,a12)') + 'Cext','Cabs','Csca','Cpr','' write(33,'(a5,5e15.6)') 'total',cext,cabs,csca,cext-cpr, + assym do i=1,nL write(33,'(i5,5e15.6)') i,cexti(i),cabsi(i),cscai(i), + cpri(i),assymi(i) enddo close(33) cz=pione*gcvr*gcvr assym=cpr/csca assymx=cprx/cscax assymy=cpry/cscay cabsxv=cabsx/cz cabsyv=cabsy/cz cextxv=cextx/cz cextyv=cexty/cz cscaxv=cscax/cz cscayv=cscay/cz cprxv=cprx/cz cprxv=cextxv-cprxv cpryv=cpry/cz cpryv=cextyv-cpryv cscav=0.5d0*(cscaxv+cscayv) cextv=0.5d0*(cextxv+cextyv) cabsv=0.5d0*(cabsxv+cabsyv) cprv=0.5d0*(cprxv+cpryv) cbakxv=cbakx/cz cbakyv=cbaky/cz cbakv=0.5d0*(cbakxv+cbakyv) temp=gcvr*gcvr/gcs cabsxs=cabsxv*temp cabsys=cabsyv*temp cextxs=cextxv*temp cextys=cextyv*temp cscaxs=cscaxv*temp cscays=cscayv*temp cprxs=cprxv*temp cprys=cpryv*temp cscas=cscav*temp cexts=cextv*temp cabss=cabsv*temp cprs=cprv*temp cbakxs=cbakxv*temp cbakys=cbakyv*temp cbaks=cbakv*temp 222 format(1x,a1,6e13.5) 221 format(6x,a5,8x,a5,8x,a5,8x,a5,8x,a4,5x,a12) write(6,221) + 'Qextv','Qabsv','Qscav','Qbakv','Qprv','' write(6,222) 't',cextv,cabsv,cscav,cbakv,cprv,assym write(6,222) 'x',cextxv,cabsxv,cscaxv,cbakxv,cprxv,assymx write(6,222) 'y',cextyv,cabsyv,cscayv,cbakyv,cpryv,assymy write(6,221) + 'Qexts','Qabss','Qscas','Qbaks','Qprs','' write(6,222) 't',cexts,cabss,cscas,cbaks,cprs,assym write(6,222) 'x',cextxs,cabsxs,cscaxs,cbakxs,cprxs,assymx write(6,222) 'y',cextys,cabsys,cscays,cbakys,cprys,assymy temp=-(cabs+csca-cext)/cext c write(6,'(/,a37,e14.5)') c + 'Accuracy of this numerical solution: ',temp open(12,file=fileout,status='unknown') write(12,'(a20,a16,a18,a4,f8.3,a5,f8.3)') fileout, + '--- input file: ',FLNAME,' xv:',xv,' xs:',xs write(12,'(a24,6f7.2)') 'Ranges of Euler angles: ', + betami,betamx,thetmi,thetmx,phaimi,phaimx write(12,'(a19,3i5,6x,a28,i5)') + 'nbeta,nthet,nphai: ',nbeta,nthet,nphai, + '# of orientations averaged: ',nram write(12,221) + 'Cext','Cabs','Csca','Cbak','Cpr','' write(12,222) + 't',cext,cabs,csca,cbak,cext-cpr,assym write(12,222) + 'x',cextx,cabsx,cscax,cbakx,cextx-cprx,assymx write(12,222) + 'y',cexty,cabsy,cscay,cbaky,cexty-cpry,assymy write(12,221) + 'Qextv','Qabsv','Qscav','Qbakv','Qprv','' write(12,222) 't',cextv,cabsv,cscav,cbakv,cprv,assym write(12,222) 'x',cextxv,cabsxv,cscaxv,cbakxv,cprxv,assymx write(12,222) 'y',cextyv,cabsyv,cscayv,cbakyv,cpryv,assymy write(12,221) + 'Qexts','Qabss','Qscas','Qbaks','Qprs','' write(12,222) 't',cexts,cabss,cscas,cbaks,cprs,assym write(12,222) 'x',cextxs,cabsxs,cscaxs,cbakxs,cprxs,assymx write(12,222) 'y',cextys,cabsys,cscays,cbakys,cprys,assymy write(12,'(/)') write(12,*) +'For the scattering plane of phi=0 (i.e., the x-z plane) only:' write(12,'(2x,a4,5x,a6,6x,a4,6x,a5,8x,a5,8x,a5,8x,a5)') + 's.a.','itotal','pol.','S1*S1','S4*S4','S3*S3','S2*S2' do i=1,nang2 write(12,'(f7.1,e13.5,f8.4,4e13.5)') + dang(i),inat(i),pol(i),i11(i),i21(i),i12(i),i22(i) enddo close(12) open(11,file='mueller.out',status='unknown') write(11,'(a11,6x,a16)') 'mueller.out','(Mueller matrix)' write(11,'(a16,a20)') 'Input filename: ',FLNAME write(11,'(a19,3i5)') 'nbeta,nthet,nphai: ', + nbeta,nthet,nphai write(11,'(a24,6f7.2)') 'Ranges of Euler angles: ', + betami,betamx,thetmi,thetmx,phaimi,phaimx write(11,'(a28,i5)') '# of orientations averaged: ',nram write(11,'(a30,f9.3)') 'interval of scattering angle: ',sang write(11,'(a30,f9.3)') 'interval of azimuth angle: ',pang do jc=1,npng t=pang*dble(jc-1) write(11,'(a18,f8.3)') 'phi (in degrees): ',t do i=1,nang2 write(11,'(f7.1,4e16.7)') + dang(i),mue(1,1,jc,i),mue(1,2,jc,i),mue(1,3,jc,i), + mue(1,4,jc,i) write(11,'(7x,4e16.7)') + mue(2,1,jc,i),mue(2,2,jc,i),mue(2,3,jc,i), + mue(2,4,jc,i) write(11,'(7x,4e16.7)') + mue(3,1,jc,i),mue(3,2,jc,i),mue(3,3,jc,i), + mue(3,4,jc,i) write(11,'(7x,4e16.7)') + mue(4,1,jc,i),mue(4,2,jc,i),mue(4,3,jc,i), + mue(4,4,jc,i) enddo enddo write(11,*) 'phi (in degrees): 360' do i=1,nang2 write(11,'(f7.1,4e16.7)') + dang(i),mue(1,1,1,i),mue(1,2,1,i),mue(1,3,1,i), + mue(1,4,1,i) write(11,'(7x,4e16.7)') + mue(2,1,1,i),mue(2,2,1,i),mue(2,3,1,i), + mue(2,4,1,i) write(11,'(7x,4e16.7)') + mue(3,1,1,i),mue(3,2,1,i),mue(3,3,1,i), + mue(3,4,1,i) write(11,'(7x,4e16.7)') + mue(4,1,1,i),mue(4,2,1,i),mue(4,3,1,i), + mue(4,4,1,i) enddo close(11) STOP END SUBROUTINE orientcd(BETAMI,BETAMX,THETMI,THETMX,PHIMIN, & PHIMAX,MXBETA,MXTHET,MXPHI,NBETA,NTHETA,NPHI, & BETA,THETA,PHI) C Arguments: INTEGER MXBETA,MXTHET,MXPHI,NBETA,NTHETA,NPHI double precision BETAMI,BETAMX,THETMI,THETMX,PHIMIN,PHIMAX double precision BETA(MXBETA),THETA(MXTHET),PHI(MXPHI) C Local variables: INTEGER J double precision DELTA C*********************************************************************** C Given: BETAMI=minimum value of beta (radians) C BETAMX=maximum value of beta (radians) C THETMI=minimum value of theta (radians) C THETMX=maximum value of theta (radians) C PHIMIN=minimum value of phi (radians) C PHIMAX=maximum value of phi (radians) C MXBETA,MXTHET,MXPHI=dimensions of the arrays BETA,THETA,PHI C NBETA=number of values of beta C NTHETA=number of values of theta C NPHI=number of values of PHI C Returns: BETA(1-NBETA)=beta values (radians) C THETA(1-NTHETA)=theta values (radians) C PHI(1-NPHI)=phi values (radians) C Purpose: to generate a sequence of desired target orientations C Present version assumes: C beta to be uniformly distributed between BETAMI and BETAMX C cos(theta) to be uniformly distributed between cos(THETMI) and C cos(THETMX) C phi to be uniformly distributed between PHIMIN and PHIMAX C If NPHI=1, first angle is THETMI, last angle is THETMX C If NPHI>1, then angles are midpoints of intervals of equal C range in theta subdividing range from THETMI to THETMX C*********************************************************************** BETA(1)=BETAMI IF(NBETA.GT.1)THEN DELTA=(BETAMX-BETAMI)/DBLE(NBETA-1) DO 1000 J=2,NBETA BETA(J)=BETA(1)+DELTA*DBLE(J-1) 1000 CONTINUE ENDIF IF(NPHI.EQ.1.AND.NTHETA.GT.1)THEN DELTA=(DCOS(THETMX)-DCOS(THETMI))/DBLE(NTHETA-1) THETA(1)=THETMI ELSE DELTA=(DCOS(THETMX)-DCOS(THETMI))/DBLE(NTHETA) THETA(1)=DACOS(DCOS(THETMI)+0.5d0*DELTA) ENDIF IF(NTHETA.GT.1)THEN DO 2000 J=2,NTHETA THETA(J)=DACOS(DCOS(THETA(1))+DELTA*DBLE(J-1)) 2000 CONTINUE ENDIF c DELTA=(PHIMAX-PHIMIN)/DBLE(NPHI) c PHI(1)=PHIMIN+0.5D0*DELTA PHI(1)=PHIMIN IF(NPHI.GT.1)THEN DELTA=(PHIMAX-PHIMIN)/DBLE(NPHI-1) DO 3000 J=2,NPHI PHI(J)=PHI(1)+DELTA*DBLE(J-1) 3000 CONTINUE ENDIF RETURN END SUBROUTINE orientud(BETAMI,BETAMX,THETMI,THETMX,PHIMIN, & PHIMAX,MXBETA,MXTHET,MXPHI,NBETA,NTHETA,NPHI, & BETA,THETA,PHI) C Arguments: INTEGER MXBETA,MXTHET,MXPHI,NBETA,NTHETA,NPHI double precision BETAMI,BETAMX,THETMI,THETMX,PHIMIN, + PHIMAX,BETA(MXBETA),THETA(MXTHET),PHI(MXPHI) C Local variables: INTEGER J double precision DELTA C*********************************************************************** C Given: BETAMI=minimum value of beta (radians) C BETAMX=maximum value of beta (radians) C THETMI=minimum value of theta (radians) C THETMX=maximum value of theta (radians) C PHIMIN=minimum value of phi (radians) C PHIMAX=maximum value of phi (radians) C MXBETA,MXTHET,MXPHI=dimensions of the arrays BETA,THETA,PHI C NBETA=number of values of beta C NTHETA=number of values of theta C NPHI=number of values of PHI C Returns: BETA(1-NBETA)=beta values (radians) C THETA(1-NTHETA)=theta values (radians) C PHI(1-NPHI)=phi values (radians) C Note: it is assumed that target orientation weight function C can be factored into WGTA*WGTB -- i.e., that rotations C around a1 are decoupled from orientation of a1. C Purpose: to generate a sequence of desired target orientations C Present version assumes: C beta to be uniformly distributed between BETAMI and BETAMX C theta to be uniformly distributed between THETMI and THETMX C phi to be uniformly distributed between PHIMIN and PHIMAX C first angle is THETMI, last angle is THETMX C*********************************************************************** BETA(1)=BETAMI IF(NBETA.GT.1)THEN DELTA=(BETAMX-BETAMI)/DBLE(NBETA-1) DO 1000 J=2,NBETA BETA(J)=BETA(1)+DELTA*DBLE(J-1) 1000 CONTINUE ENDIF THETA(1)=THETMI IF(NTHETA.GT.1)THEN DELTA=(THETMX-THETMI)/DBLE(NTHETA-1) DO 2000 J=2,NTHETA THETA(J)=THETA(1)+DELTA*DBLE(J-1) 2000 CONTINUE ENDIF PHI(1)=PHIMIN IF(NPHI.GT.1)THEN DELTA=(PHIMAX-PHIMIN)/DBLE(NPHI-1) DO 3000 J=2,NPHI PHI(J)=PHI(1)+DELTA*DBLE(J-1) 3000 CONTINUE ENDIF RETURN END DOUBLE PRECISION FUNCTION ran1d(idum) INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV DOUBLE PRECISION AM,EPS,RNMX PARAMETER (IA=16807,IM=2147483647,AM=1.d0/IM, *IQ=127773,IR=2836,NTAB=32,NDIV=1+(IM-1)/NTAB, *EPS=1.2d-7,RNMX=1.d0-EPS) INTEGER j,k,iv(NTAB),iy SAVE iv,iy DATA iv /NTAB*0/, iy /0/ if (idum.le.0.or.iy.eq.0) then idum=max(-idum,1) do 11 j=NTAB+8,1,-1 k=idum/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM if (j.le.NTAB) iv(j)=idum 11 continue iy=iv(1) endif k=idum/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM j=1+iy/NDIV iy=iv(j) iv(j)=idum ran1d=dmin1(AM*iy,RNMX) 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 common/cofmnv0/cof0(ni0) common/crot/cofsr(nmp) common/fnr/fnr(0:2*(np+2)) 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:2*(np+2)) 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) 1 *fnr(n+k-1)*dkt)/(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:2*(np+2)) 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 subroutine mueller.f c returns the values of 16 elements of the 4x4 c Mueller matrix calculated from the known 2x2 c amplitude scattering matrix c using Bohren & Huffman's formulas (p.65) subroutine mueller(s1,s2,s3,s4,s) double precision s(4,4),s1s,s2s,s3s,s4s complex*16 s1,s2,s3,s4 complex*16 s2s3c,s1s4c,s2s4c,s1s3c,s1s2c complex*16 s3s4c,s2s1c,s4s3c,s2cs4,s3cs1 s1s=cdabs(s1)**2 s2s=cdabs(s2)**2 s3s=cdabs(s3)**2 s4s=cdabs(s4)**2 s2s3c=s2*dconjg(s3) s1s4c=s1*dconjg(s4) s2s4c=s2*dconjg(s4) s1s3c=s1*dconjg(s3) s1s2c=s1*dconjg(s2) s3s4c=s3*dconjg(s4) s2s1c=s2*dconjg(s1) s4s3c=s4*dconjg(s3) s2cs4=dconjg(s2)*s4 s3cs1=dconjg(s3)*s1 s(1,1)=0.5d0*(s1s+s2s+s3s+s4s) s(1,2)=0.5d0*(s2s-s1s+s4s-s3s) s(1,3)=s2s3c+s1s4c s(1,4)=dimag(s2s3c-s1s4c) s(2,1)=0.5d0*(s2s-s1s-s4s+s3s) s(2,2)=0.5d0*(s2s+s1s-s4s-s3s) s(2,3)=s2s3c-s1s4c s(2,4)=dimag(s2s3c+s1s4c) s(3,1)=s2s4c+s1s3c s(3,2)=s2s4c-s1s3c s(3,3)=s1s2c+s3s4c s(3,4)=dimag(s2s1c+s4s3c) s(4,1)=dimag(s2cs4+s3cs1) s(4,2)=dimag(s2cs4-s3cs1) s(4,3)=dimag(s1s2c-s3s4c) s(4,4)=s1s2c-s3s4c 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 subroutine rtr(anpt,nodrj,nodri,ekt,drot) 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) double precision drot(nrc) complex*16 anpt(2,nmp),ant(2,2*np),amt(2,-np:np),a,b, + ekt(np),eks(-np:np),atr(2,np,nmp) common/tran/atr c eks(0)=1.d0 nmax=max(nodrj,nodri) do m=1,nmax eks(m)=ekt(m) eks(-m)=dconjg(eks(m)) enddo irc=0 do n=1,nodrj n1=n*(n+1) do m=-n,n amt(1,m)=0.d0 amt(2,m)=0.d0 enddo do k=-n,n kn=n1+k a=eks(k)*anpt(1,kn) b=eks(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 enddo do m=-n,n imn=n1+m anpt(1,imn)=amt(1,m) anpt(2,imn)=amt(2,m) enddo enddo mmax=min(nodrj,nodri) do m=-mmax,mmax n1=max(1,iabs(m)) do n=n1,nodrj imn=n*(n+1)+m do ip=1,2 ant(ip,n)=anpt(ip,imn) enddo enddo do n=n1,nodri imn=n*(n+1)+m a=0.d0 b=0.d0 do l=n1,nodrj 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) enddo anpt(1,imn) = a anpt(2,imn) = b enddo enddo in=1 irc=0 do n=1,nodri in=-in n1=n*(n+1) do m=-n,n amt(1,m)=0.d0 amt(2,m)=0.d0 enddo sik=-in do 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 enddo sik=-in do m=-n,n sik=-sik imn=n1+m anpt(1,imn)=amt(1,m)*eks(-m)*sik anpt(2,imn)=amt(2,m)*eks(-m)*sik enddo enddo return end subroutine trv(anpt,nodrj,nodri) 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 mmax=min(nodrj,nodri) do m=-mmax,mmax n1=max(1,iabs(m)) do n=n1,nodrj imn=n*(n+1)+m do ip=1,2 ant(ip,n)=anpt(ip,imn) enddo enddo do n=n1,nodri imn=n*(n+1)+m a=0.d0 b=0.d0 do l=n1,nodrj 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) enddo anpt(1,imn) = a anpt(2,imn) = b enddo enddo 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 subroutine transs(k,nL,r0,nmax0,nmax,uvmax,fint,ek, + besj,besy,drot,as,bs,as1,bs1,ind,icase) implicit double precision (a-h,o-z) include 'gmm01f.par' parameter (nmp=np*(np+2),nmp0=(np+1)*(np+4)/2) parameter (ni0=np*(np+1)*(2*np+1)/3+np*np) parameter (nrc=4*np*(np+1)*(np+2)/3+np) parameter (ng0=np*(2*np**3+10*np**2+19*np+5)/6) integer u,v,nmax(nLp),uvmax(nLp),ind(nLp) double precision r0(9,nLp),drot(nrc),k, + besj(0:2*np+1),besy(0:2*np+1) complex*16 atr(2,np,nmp), + ek(np),as(nLp,nmp),bs(nLp,nmp),as1(nLp,nmp), + bs1(nLp,nmp),at1(2,nmp),t1,t2,ephi common/tran/atr common/rot/bcof(0:np+2),dc(-np:np,0:nmp) common/fnr/fnr(0:2*(np+2)) common/pitau/pi(nmp0),tau(nmp0) common/ig0/iga0(ni0) common/g0/ga0(ng0) common/cofmnv0/cof0(ni0) common/crot/cofsr(nmp) 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 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.le.1.d0) goto 1 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 file' write(6,*) + 'for the aggregate configuration and then' write(6,*) 'run again.' stop endif 1 if(temp.le.fint) goto 10 ephi=dcmplx(cphi,sphi) nlarge=max(nmax(i),nmax(j)) do m=1,nlarge ek(m)=ephi**m enddo xd=k*d nbes=2*nlarge+1 call besseljd(nbes,xd,besj) call besselyd(nbes,xd,besy) 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)=dc(u,imn) enddo enddo enddo nsmall=min(nmax(i),nmax(j)) 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) do m=-nsmall,nsmall n1=max(1,iabs(m)) do n=n1,nlarge do v=n1,nlarge iuv=v*(v+1)+m if(icase.lt.2) then call cofxuds0(nmax0,m,n,v,besj,besy, + atr(1,n,iuv),atr(2,n,iuv),t1,t2) else call cofxuds0(nmax0,m,n,v,besj,besy, + t1,t2,atr(1,n,iuv),atr(2,n,iuv)) endif if(x0.eq.0.d0.and.y0.eq.0.d0) then if(z0.lt.0.d0) then sic=dble((-1)**(n+v)) atr(1,n,iuv)=sic*atr(1,n,iuv) atr(2,n,iuv)=-sic*atr(2,n,iuv) endif endif 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 if(x0.eq.0.d0.and.y0.eq.0.d0) then call trv(at1,nmax(j),nmax(i)) else call rtr(at1,nmax(j),nmax(i),ek,drot) endif do imn=1,uvmax(i) 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 if(x0.eq.0.d0.and.y0.eq.0.d0) then call trv(at1,nmax(i),nmax(j)) else call rtr(at1,nmax(i),nmax(j),ek,drot) endif do imn=1,uvmax(j) 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 C C subroutine scoatabd.f C slightly edited from SCSERD.FOR BY R.T. Wang C SUB SCSERD: DP CORE-MANTLE-SPHERE MIE COEFFICIENTS C 6-7-91 Ru T. Wang C ED. OF SCSE1W(1974(UNIVAC) & 10-84(VAX)) & W.X.WANG(5-7-85). C DOUBLE PRECISION,DOWNWARD RECURSION RATIO FUNCTION ALGORITHMS. C Q=(CORE RADIUS)/(MANTLE RADIUS); XB=SIZE PARAMETER OF MANTLE C COMPLEX REFRACTIVE INDEX OF CORE = (XM1,YM1) C COMPLEX REFRACTIVE INDEX OF MANTLE = (XM2,YM2) C SUBROUTINE scoatabd(XB,Q,XM1,YM1,XM2,YM2,np,an,bn, + NADD,NSTOP) IMPLICIT double precision (A-H,P-Z) parameter (nab=500,ndx=5000) dimension AR(nab),AI(nab),BR(nab),BI(nab), 1 AM1AR(ndx),AM1AI(ndx),AM2AR(ndx),AM2AI(ndx), 2 AM2BR(ndx),AM2BI(ndx),AB(ndx),SM2AR(ndx), 3 SM2AI(ndx),SM2BR(ndx),SM2BI(ndx),SB(ndx), 4 BM2AR(ndx),BM2AI(ndx),BM2BR(ndx),BM2BI(ndx), 5 BDBR(ndx),BDBI(ndx),BB(ndx),CM2AR(ndx), 6 CM2AI(ndx),CM2BR(ndx),CM2BI(ndx),CB(ndx), 7 U(4),V(4),D1(4),EV(4),SHV(4),CHV(4),SU(4), 8 CU(4),SN1R(4),SN1I(4),CN1R(4),CN1I(4) complex*16 an(np),bn(np) C do i=1,np an(i)=0.d0 bn(i)=0.d0 enddo XM=DMAX1(XM1,XM2) YM1P=DABS(YM1) YM2P=DABS(YM2) YM=DMAX1(YM1P,YM2P) XN=XB*DSQRT(XM**2+YM**2) NX=1.1D0*XN+10.D0 if(NX.gt.ndx) then write(6,*) + 'parameter (ndx) in sub. scoatabd too small' write(6,*) 'please change ndx to ',NX write(6,*) ' recompile, and then try again' stop endif oneth=1.d0/3.d0 NSTOP=XB+4.d0*XB**oneth NSTOP=NSTOP+2+NADD if(NSTOP.gt.ndx) then write(6,*) 'particle size too large' stop endif XA=XB*Q CTST=1.0D-14 U(1)=XM1*XA V(1)=YM1*XA U(2)=XM2*XA V(2)=YM2*XA U(3)=XM2*XB V(3)=YM2*XB U(4)=XB V(4)=0.0D0 K=1 IF(Q.EQ.0.0D0) K=3 DO 10 J=K,4 D1(J)=U(J)*U(J)+V(J)*V(J) EV(J)=DEXP(V(J)) SHV(J)=0.5D0*(EV(J)*EV(J)-1.0D0) CHV(J)=0.5D0*(EV(J)*EV(J)+1.0D0) SU(J)=DSIN(U(J)) CU(J)=DCOS(U(J)) SN1R(J)=SU(J)*CHV(J) SN1I(J)=CU(J)*SHV(J) CN1R(J)=CU(J)*CHV(J) 10 CN1I(J)=-SU(J)*SHV(J) CNX=DFLOAT(NX) FCT0=2.0D0*CNX+3.0D0 IF(Q.EQ.0.0D0) GO TO 12 SM1A0R=U(1)/FCT0 SM1A0I=V(1)/FCT0 SM2A0R=U(2)/FCT0 SM2A0I=V(2)/FCT0 12 SM2B0R=U(3)/FCT0 SM2B0I=V(3)/FCT0 SB0=U(4)/FCT0 DO 15 I=1,NX N=NX-I+1 CNN=DFLOAT(N) FCT=2.0D0*CNN+1.0D0 IF(Q.EQ.0.0D0) GO TO 18 SM1A1R=+U(1)*FCT/D1(1)-SM1A0R SM1A1I=-V(1)*FCT/D1(1)-SM1A0I SM2A1R=+U(2)*FCT/D1(2)-SM2A0R SM2A1I=-V(2)*FCT/D1(2)-SM2A0I AM1AR(N)=-CNN*U(1)/D1(1)+SM1A1R AM1AI(N)=+CNN*V(1)/D1(1)+SM1A1I AM2AR(N)=-CNN*U(2)/D1(2)+SM2A1R AM2AI(N)=+CNN*V(2)/D1(2)+SM2A1I QSM1A=SM1A1R**2+SM1A1I**2 QSM2A=SM2A1R**2+SM2A1I**2 SM1A0R=+SM1A1R/QSM1A SM1A0I=-SM1A1I/QSM1A SM2A0R=+SM2A1R/QSM2A SM2A0I=-SM2A1I/QSM2A SM2AR(N)=SM2A0R SM2AI(N)=SM2A0I 18 SM2BR(N)=+U(3)*FCT/D1(3)-SM2B0R SM2BI(N)=-V(3)*FCT/D1(3)-SM2B0I AM2BR(N)=-CNN*U(3)/D1(3)+SM2BR(N) AM2BI(N)=+CNN*V(3)/D1(3)+SM2BI(N) QSM2B=SM2BR(N)**2+SM2BI(N)**2 SM2B0R=+SM2BR(N)/QSM2B SM2B0I=-SM2BI(N)/QSM2B AB(N)=(CNN+1.0D0)/U(4)-SB0 SB0=U(4)/(FCT-SB0*U(4)) 15 SB(N)=SB0 IF(Q.EQ.0.0D0) GO TO 20 QCM2A0=CN1R(2)**2+CN1I(2)**2 BM2A0R=-SN1R(2)*CN1R(2)-SN1I(2)*CN1I(2) BM2A0R=BM2A0R/QCM2A0 BM2A0I=+SN1R(2)*CN1I(2)-SN1I(2)*CN1R(2) BM2A0I=BM2A0I/QCM2A0 20 QSM2B0=SN1R(3)**2+SN1I(3)**2 QCM2B0=CN1R(3)**2+CN1I(3)**2 BM2B0R=-SN1R(3)*CN1R(3)-SN1I(3)*CN1I(3) BM2B0R=BM2B0R/QCM2B0 BM2B0I=+SN1R(3)*CN1I(3)-SN1I(3)*CN1R(3) BM2B0I=BM2B0I/QCM2B0 BB0=-SN1R(4)/CN1R(4) BDB0R=0.0D0 BDB0I=-1.0D0 IF(Q.EQ.0.0D0) GO TO 22 UM2R0=+SN1R(2)*SN1R(3)+SN1I(2)*SN1I(3) UM2R0=UM2R0/QSM2B0 UM2I0=-SN1R(2)*SN1I(3)+SN1I(2)*SN1R(3) UM2I0=UM2I0/QSM2B0 VM2R0=+CN1R(2)*CN1R(3)+CN1I(2)*CN1I(3) VM2R0=VM2R0/QCM2A0 VM2I0=+CN1R(2)*CN1I(3)-CN1I(2)*CN1R(3) VM2I0=VM2I0/QCM2A0 22 CONTINUE DO 25 N=1,NX CN=DFLOAT(N) IF(Q.EQ.0.0D0) GO TO 24 CM2A0R=+CN*U(2)/D1(2)-BM2A0R CM2A0I=-CN*V(2)/D1(2)-BM2A0I QCM2A=CM2A0R**2+CM2A0I**2 CM2AR(N)=+CM2A0R/QCM2A CM2AI(N)=-CM2A0I/QCM2A BM2AR(N)=-CN*U(2)/D1(2)+CM2AR(N) BM2AI(N)=+CN*V(2)/D1(2)+CM2AI(N) BM2A0R=BM2AR(N) BM2A0I=BM2AI(N) 24 CM2BR(N)=+CN*U(3)/D1(3)-BM2B0R CM2BI(N)=-CN*V(3)/D1(3)-BM2B0I CB(N)=CN/U(4)-BB0 DBDB0R=+CN/U(4)-BDB0R DBDB0I=-BDB0I QCM2B=CM2BR(N)**2+CM2BI(N)**2 QDB=DBDB0R**2+DBDB0I**2 CM2B1R=+CM2BR(N)/QCM2B CM2B1I=-CM2BI(N)/QCM2B CB1=1.0D0/CB(N) DBDB1R=DBDB0R/QDB DBDB1I=-DBDB0I/QDB BM2BR(N)=-CN*U(3)/D1(3)+CM2B1R BM2BI(N)=+CN*V(3)/D1(3)+CM2B1I BM2B0R=BM2BR(N) BM2B0I=BM2BI(N) BB(N)=-CN/U(4)+CB1 BB0=BB(N) BDBR(N)=-CN/U(4)+DBDB1R BDBI(N)=DBDB1I BDB0R=BDBR(N) BDB0I=BDBI(N) SNB=SB(N)*SN1R(4) CNB=CB(N)*CN1R(4) SSB1R=AM2BR(N)-XM2*AB(N) SSB1I=AM2BI(N)-YM2*AB(N) SZB1R=AM2BR(N)-XM2*BDBR(N)+YM2*BDBI(N) SZB1I=AM2BI(N)-YM2*BDBR(N)-XM2*BDBI(N) SSB2R=XM2*AM2BR(N)-YM2*AM2BI(N)-AB(N) SSB2I=XM2*AM2BI(N)+YM2*AM2BR(N) SZB2R=XM2*AM2BR(N)-YM2*AM2BI(N)-BDBR(N) SZB2I=XM2*AM2BI(N)+YM2*AM2BR(N)-BDBI(N) IF(Q.NE.0.0D0) GO TO 45 ANNR=SNB*SSB1R ANNI=SNB*SSB1I ANDR=SNB*SZB1R-CNB*SZB1I ANDI=CNB*SZB1R+SNB*SZB1I BNNR=SNB*SSB2R BNNI=SNB*SSB2I BNDR=SNB*SZB2R-CNB*SZB2I BNDI=CNB*SZB2R+SNB*SZB2I GO TO 65 45 UM2R1=SM2AR(N)*SM2BR(N)-SM2AI(N)*SM2BI(N) UM2I1=SM2AR(N)*SM2BI(N)+SM2AI(N)*SM2BR(N) UM2R=UM2R1*UM2R0-UM2I1*UM2I0 UM2I=UM2R1*UM2I0+UM2I1*UM2R0 VM2R1=CM2AR(N)*CM2BR(N)-CM2AI(N)*CM2BI(N) VM2I1=CM2AR(N)*CM2BI(N)+CM2AI(N)*CM2BR(N) VM2R=VM2R1*VM2R0-VM2I1*VM2I0 VM2I=VM2R1*VM2I0+VM2I1*VM2R0 UVR=UM2R*VM2R-UM2I*VM2I UVI=UM2R*VM2I+UM2I*VM2R SSA1R=XM1*AM2AR(N)-YM1*AM2AI(N) 1 -XM2*AM1AR(N)+YM2*AM1AI(N) SSA1I=YM1*AM2AR(N)+XM1*AM2AI(N) 1 -YM2*AM1AR(N)-XM2*AM1AI(N) CSB1R=BM2BR(N)-XM2*AB(N) CSB1I=BM2BI(N)-YM2*AB(N) SCS1R=SSA1R*CSB1R-SSA1I*CSB1I SCS1I=SSA1R*CSB1I+SSA1I*CSB1R UVS1R=UVR*SCS1R-UVI*SCS1I UVS1I=UVR*SCS1I+UVI*SCS1R CSA1R=XM1*BM2AR(N)-YM1*BM2AI(N) 1 -XM2*AM1AR(N)+YM2*AM1AI(N) CSA1I=YM1*BM2AR(N)+XM1*BM2AI(N) 1 -YM2*AM1AR(N)-XM2*AM1AI(N) CSS1R=CSA1R*SSB1R-CSA1I*SSB1I CSS1I=CSA1R*SSB1I+CSA1I*SSB1R CZB1R=BM2BR(N)-XM2*BDBR(N)+YM2*BDBI(N) CZB1I=BM2BI(N)-YM2*BDBR(N)-XM2*BDBI(N) SCZ1R=SSA1R*CZB1R-SSA1I*CZB1I SCZ1I=SSA1R*CZB1I+SSA1I*CZB1R UVZ1R=UVR*SCZ1R-UVI*SCZ1I UVZ1I=UVI*SCZ1R+UVR*SCZ1I CSZ1R=CSA1R*SZB1R-CSA1I*SZB1I CSZ1I=CSA1R*SZB1I+CSA1I*SZB1R SSA2R=XM2*AM2AR(N)-YM2*AM2AI(N) 1 -XM1*AM1AR(N)+YM1*AM1AI(N) SSA2I=YM2*AM2AR(N)+XM2*AM2AI(N) 1 -YM1*AM1AR(N)-XM1*AM1AI(N) CSB2R=XM2*BM2BR(N)-YM2*BM2BI(N)-AB(N) CSB2I=XM2*BM2BI(N)+YM2*BM2BR(N) SCS2R=SSA2R*CSB2R-SSA2I*CSB2I SCS2I=SSA2R*CSB2I+SSA2I*CSB2R UVS2R=UVR*SCS2R-UVI*SCS2I UVS2I=UVR*SCS2I+UVI*SCS2R CSA2R=XM2*BM2AR(N)-YM2*BM2AI(N) 1 -XM1*AM1AR(N)+YM1*AM1AI(N) CSA2I=YM2*BM2AR(N)+XM2*BM2AI(N) 1 -YM1*AM1AR(N)-XM1*AM1AI(N) CSS2R=CSA2R*SSB2R-CSA2I*SSB2I CSS2I=CSA2R*SSB2I+CSA2I*SSB2R CZB2R=XM2*BM2BR(N)-YM2*BM2BI(N)-BDBR(N) CZB2I=XM2*BM2BI(N)+YM2*BM2BR(N)-BDBI(N) SCZ2R=SSA2R*CZB2R-SSA2I*CZB2I SCZ2I=SSA2R*CZB2I+SSA2I*CZB2R UVZ2R=UVR*SCZ2R-UVI*SCZ2I UVZ2I=UVI*SCZ2R+UVR*SCZ2I CSZ2R=CSA2R*SZB2R-CSA2I*SZB2I CSZ2I=CSA2R*SZB2I+CSA2I*SZB2R 55 ANNR=SNB*(UVS1R-CSS1R) ANNI=SNB*(UVS1I-CSS1I) ANDR=SNB*(UVZ1R-CSZ1R)-CNB*(UVZ1I-CSZ1I) ANDI=CNB*(UVZ1R-CSZ1R)+SNB*(UVZ1I-CSZ1I) BNNR=SNB*(UVS2R-CSS2R) BNNI=SNB*(UVS2I-CSS2I) BNDR=SNB*(UVZ2R-CSZ2R)-CNB*(UVZ2I-CSZ2I) BNDI=CNB*(UVZ2R-CSZ2R)+SNB*(UVZ2I-CSZ2I) 65 AND=ANDR*ANDR+ANDI*ANDI IF(AND.NE.0.0D0) GO TO 70 ABANDR=DABS(ANDR) ABANDI=DABS(ANDI) ABANNR=DABS(ANNR) ABANNI=DABS(ANNI) AAAA=DMAX1(ABANDR,ABANDI,ABANNR,ABANNI) ANDR=ANDR/AAAA ANDI=ANDI/AAAA ANNR=ANNR/AAAA ANNI=ANNI/AAAA GO TO 65 70 CONTINUE AR(N)=(ANNR*ANDR+ANNI*ANDI)/AND AI(N)=(ANNI*ANDR-ANNR*ANDI)/AND 75 BND=BNDR*BNDR+BNDI*BNDI IF(BND.NE.0.0D0) GO TO 80 ABBNDR=DABS(BNDR) ABBNDI=DABS(BNDI) ABBNNR=DABS(BNNR) ABBNNI=DABS(BNNI) BBBB=DMAX1(ABBNDR,ABBNDI,ABBNNR,ABBNNI) BNDR=BNDR/BBBB BNDI=BNDI/BBBB BNNR=BNNR/BBBB BNNI=BNNI/BBBB GO TO 75 80 CONTINUE BR(N)=(BNNR*BNDR+BNNI*BNDI)/BND BI(N)=(BNNI*BNDR-BNNR*BNDI)/BND TI=AR(N)*AR(N)+AI(N)*AI(N)+BR(N)*BR(N) 1 +BI(N)*BI(N) TI=TI/(AR(1)*AR(1)+AI(1)*AI(1)+BR(1)*BR(1)+ 1 BI(1)*BI(1)) SN1R(4)=SNB CN1R(4)=CNB UM2R0=UM2R UM2I0=UM2I VM2R0=VM2R VM2I0=VM2I IF(NSTOP-N) 50,50,25 C IF(TI-CTST) 50,50,25 25 CONTINUE 50 CONTINUE do i=1,NSTOP an(i)=dcmplx(AR(i),-AI(i)) bn(i)=dcmplx(BR(i),-BI(i)) enddo RETURN END