C C DOUBLE PRECISION Fortran code gmm03fq.f for calculation of elastic C radiative scattering by an ensemble of variously shaped small C particles in a single fixed-orientation C In all aspects and usage this gmm03fq.f code is exactly the same as C the gmm03f.f code except that the double precision subroutine tm0d.f C in the gmm03f.f code is replaced by the extended precision subroutine C tm0q.f in this code C Yu-lin Xu C Released to public 6/2003 C available online at http://www.astro.ufl.edu/~xu 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 (1) The multiparticle-scattering formulation used in this code can be C found in the papers C Bruning and Lo, IEEE Trans. Anten. Prop. AP-19, 378 (1971) C Xu, Appl. Opt. 34, 4573 (1995) C Appl. Opt. 36, 9496 (1997) C Phys. Lett. A 249, 30 (1998) C Phys. Rev. E 67, 046620 (2003) C J. Opt. Soc. Am. A 20, 2093 (2003) C Xu and Wang, Phys. Rev. E 58, 3931 (1998) C Xu, Gustafson, Giovane, Blum, and Tehranian, C Phys. Rev. E 60, 2347 (1999) C Xu and Gustafson, JQSRT 70, 395 (2001) C Xu and Khlebtsov, JQSRT 79-80, 1121 (2003) C (2) Numerical techniques used in this code can be found in the papers C Stein, Q. Appl. Math. 19, 15 (1961) C Cruzan, Q. Appl. Math. 20, 33 (1962) C Mackowski, Proc. R. Soc. Lond. A 433, 599 (1991) C Wang and van de Hulst, Appl. Opt. 30, 106 (1991) C H.A. van der Vorst, SIAM J. Sci. Stat. Comput. 13, 631 (1992) C Gutknecht, SIAM J. Sci. Comput. 14, 1020 (1993) C Xu, J. Comput. Appl. Math. 85, 53 (1997) C J. Comput. Phys. 139, 137 (1998) C (3) In the input of this code individual particles can be an C arbitrary mixture of homogeneous spheres, core-mantle spheres, C and certain types of axially symmetric shapes including spheroids, C finite circular cylinders, Chebyshev particles, and generalized C Chebyshev particles (simulating the shapes of distorted water C drops). The subroutine tm0q.f together with its auxiliaries in C this gmm03fq.f code are a part of Mishchenko's public-domain code C tmq.new.f available at www.giss.nasa.gov/~crmim. It computes the C T-matrix of an individual nonspherical particle with an axially C symmetric shape in single-body scattering, based on Waterman's C extended boundary condition method (or called the null field C method). C C C************************* A SPECIAL NOTE ***************************C C Prior to running this code, a user must check carefully the C C configuration of an ensemble of particles to be calculated and C C make sure that all component particles are not overlapping. C C Unlike this author's other codes such as gmm01f.f, gmm01s.f, C C gmm02f.f, gmm02s.f, gmm01TrA.f, and gmm02TrA.f that are for C C ensembles of spherical particles, this gmm03fq.f code does not C C inspect the input data file for the configuration of an ensemble C C of particles. An incorrect configuration (that has some component C C particles intersecting) will result in misleading numerical C C solutions or cause divergence. It is users' responsibility to C C make sure that the input for the configuration of an ensemble of C C prticles is valid, i.e., no component particles are intersecting. C C Users must also make sure that the T-matrices of all nonspherical C C component particles in single-body scattering can be calculated by C C the subroutine tm0q.f (i.e., Mishchenko's original code tmq.new.f). C C**********************************************************************C C C C Sample parameter and input files for this code C*********************************************************************** C This example is for a single, prolate finite circular cylinder. Its C axis of rotational symmetry is aligned along the x axis. C C (1) parameter file of gmm01f.par C ---------------------------- C parameter (nLp=10,np=20) C C (2) parameter file of tmq.par.f C --------------------------- C PARAMETER (NPN1=140, NPNG1=500, NPNG2=2*NPNG1, NPN2=2*NPN1, C + NPL=NPN2+1, NPN3=NPN1+1, C + NPN4=100, NPN5=2*NPN4, NPN6=NPN4+1, NPL1=NPN5+1) C C (3) two input files C --------------- C (gmm03s.in) C gmm03.cfg C 0 C 0 0 888 C 1.d-8 600 C 0 C 0. C 1. 1. C C (gmm03.cfg) C 3.1416 1 C 1 C -2 0.5 0. 90. 0. 0. 0. 3. 1.61 0.004 0. 0. 0. C C (4) There are six output files: C gmm03fq.out (cross sections, efficiencies, and S1*S1, S2*S2, C S3*S3, S4*S4 for the scattering plane of phi=0, C i.e., the x-z plane, where S1,S2,S3,S4 are the C amplitude scattering matrix elements) C crgmm03fq.out (total and individual-particle cross sections) C amp.forw.out (the amplitude matrix for the front hemisphere) C amp.back.out (the amplitude matrix for the back hemisphere) C mue.forw.out (Mueller matrix for the front hemisphere) C mue.back.out (Mueller matrix for the back hemisphere) C C C----------------------------------------------------------------------- PROGRAM gmm03fq implicit double precision (a-h,o-z) include 'gmm01f.par' include 'tmq.par.f' C----------------------------------------------------------------------- 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": C parameter (nLp=100,np=20) C C The parameter file of "tmq.par.f" is required by the subroutine C tm0q.f. Please refer to Mishchenko's original tmq.new.f code for C details. An example for "tmq.par.f": C PARAMETER (NPN1=140, NPNG1=500, NPNG2=2*NPNG1, NPN2=2*NPN1, C + NPL=NPN2+1, NPN3=NPN1+1, C + NPN4=100, NPN5=2*NPN4, NPN6=NPN4+1, NPL1=NPN5+1) C C C C ********** Computer memory required by this code is at the level of C nLp*np^4 C This gmm03fq.f code places no restriction on the overall dimension of C an ensemble. For spherical component particles, an individual sphere C may reach a size parameter way beyond ~200. Users need to check the C original "ampld.new.f" code by Mishchenko for restrictions on size, C refractive index (especially the imaginary part), and aspect ratio of C nonspherical component particles. The total number of particles in an C ensemble that can be handled depends on the available amount of C computer memory and the largest size and/or aspect ratio of the C individual component particles. C C C----------------------------------------------------------------------- parameter (nangmax=1801) parameter (nmp=np*(np+2),nmp0=(np+1)*(np+4)/2) 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),idshp(nLp) double precision k,r0(9,nLp),x(nLp),xc(nLp),r00(3,nLp), + besj(0:2*np+1),besy(0:2*np+1),smue(4,4),i11(nangmax), + i21(nangmax),i22(nangmax),i12(nangmax),inat(nangmax), + pol(nangmax),dang(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) double precision shp(3,nLp), & RT11(NPN6,NPN4,NPN4),RT12(NPN6,NPN4,NPN4), & RT21(NPN6,NPN4,NPN4),RT22(NPN6,NPN4,NPN4), & IT11(NPN6,NPN4,NPN4),IT12(NPN6,NPN4,NPN4), & IT21(NPN6,NPN4,NPN4),IT22(NPN6,NPN4,NPN4) complex*16 A,B,cmz,Aj,Bj,A2,B2,Aj2,Bj2,A0,B0,ci,cin, + s2x,s4x,s3y,s1y,s2xj,s4xj,s3yj,s1yj, + ref(nLp),refc(nLp),atr(2,np,nmp),at(nmp),bt(nmp),ek(np), + p0(nLp,nmp),q0(nLp,nmp),an(np),bn(np),atj(nmp),btj(nmp), + as00(nLp,nmp),bs00(nLp,nmp),as0(nLp,nmp),bs0(nLp,nmp), + as(nLp,nmp),bs(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), + asx(nLp,nmp),bsx(nLp,nmp),aty(nmp),bty(nmp),atyj(nmp), + btyj(nmp),Ay(2),By(2),Ayj(2),Byj(2), + tbar(nLp,2,2,nmp,nmp),tbar0(2,2,nmp,nmp),bar(2,2), + ekt(-2*np:2*np) CHARACTER FLNAME*20,fileout*20,fileout1*19,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) COMMON /TMAT/ RT11,RT12,RT21,RT22,IT11,IT12,IT21,IT22 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 eps=1.d-20 OPEN(UNIT=1,FILE='gmm03s.in',STATUS='OLD') READ(1,'(a20)') FLNAME C----------------------------------------------------------------------- C FLNAME: input file name for the configuration of an ensemble C the first line in this file is the incident wavelength, C the second line is the total number of particles, C rest lines specify shape, orientation, position of C particle-center, size, and complex refractive index for C each component particles; one line for a particle and C each line must contain thirteen numbers, including C idshp, shp(1), shp(2), shp(3), C x, y, z, r, Re(m), Im(m), rc, Rec(m), Imc(m) C idshp -- an integer number to specify particle shape C (the same as NP in Mishchenko's code) C 0 for a sphere C -1 for a spheroid C -2 for a finite cylinder C -3 for generalized Chebyshev particles (describing C the shape of distorted water drops) C idshp>0 for Chebyshev particles, which is the degree C of the Chebyshev polynomial C shp(1) -- to specify the aspect ratio of a nonspherical C particle (the parameter EPS in Mishchenko's C "ampld.new.f" code), C for a spheroid, it is the ratio of the horizontal to C rotational axises C for a finite circular cylinder, it is the C diameter-to-length ratio C for Chebyshev particles, it is the deformation C parameter C for spheres, it is ineffective and can be set to 0 C shp(2) and shp(3) -- the azimuth and zenith angles of the C axis of rotational symmetry of a nonspherical particle C to specify its orientation, i.e., the Euler angles of C alpha and beta to rotate the coordinate system so that C the z axis after rotation is along the axis of C rotational symmetry of the particle C for spheres, simply put shp(2)=shp(3)=0 C x,y,z -- the Cartesian coordinates of a particle-center C r -- equivalent-sphere radius to specify particle size, C for spherical particles, it is the real sphere radius C Re(m), Im(m) -- the real and imaginary parts of the C refractive index of a particle C (for a core-mantle sphere, these refer to the mantle) C rc, Rec(m), Imc(m) -- effective for a core-mantle sphere C only, which are the radius, real and imaginary parts C of the refractive index of the core of a core-mantle C sphere; for all other types of particles these three C parameters are ineffective and can simply be set to 0 C C----------------------------------------------------------------------- write(6,'(a12,a20)') 'input file: ',FLNAME READ(1,*) idMie C----------------------------------------------------------------------- C idMie=1: calculating only coherent scattering, no interaction C----------------------------------------------------------------------- write(6,'(a7,i3)') 'idMie: ',idMie READ(1,*) idalpha,idbeta,iseed C----------------------------------------------------------------------- C In the input file for an ensemble configuration, a user must specify C the orientation of a nonspherical component particle, i.e., shp(2) C and shp(3), the azimuth and zenith angles of the axis of rotational C symmetry. When idalpha=1, the user-specified shp(2) will be replaced C by values randomly picked up by this code. Similarly, when idbeta=1, C the input shp(3) will be replaced by random values. Otherwise, the C user-specified values of shp(2) and shp(3) in the input will be used. C The quantity of "iseed" is an arbitrary integer, which is the seed C number used in the random number generator, effective when idalpha=1 C or idbeta=1. C----------------------------------------------------------------------- write(6,'(a,2i4,i8)') 'idalpha,idbeta,iseed: ', + idalpha,idbeta,iseed if(idalpha.eq.1) then write(6,*) + 'Note that the azimuth angle of the axis of rotational' write(6,*) + 'symmetry of nonspherical component particles are' write(6,*) + 'randomly picked up by this code' endif if(idbeta.eq.1) then write(6,*) + 'Note that the zenith angle of the axis of rotational' write(6,*) + 'symmetry of nonspherical component particles are' write(6,*) + 'randomly picked up by this code' endif READ(1,*) small,MXINT C----------------------------------------------------------------------- C small: error tolerance in the iterative process of solving boundary C condition equations (e.g., 1.d-8) C MXINT is the maximum number of iterations allowed in the iterative C process of solving boundary equations C This code uses BI-CGSTAB, the stabilized Bi-Conjugate Gradient [see C H.A. van der Vorst, SIAM J. Sci. Stat. Comput. 13, 631, (1992); C Gutknecht, SIAM J. Sci. Comput. 14, 1020 (1993)] C----------------------------------------------------------------------- write(6,'(a22,e10.2)') 'Convergence criterion:',small 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 Note that NADD is effective for spherical component particles only C Normally, set NADD=0 C----------------------------------------------------------------------- write(6,'(a,i4)') + 'Scattering orders added to Wiscombe criterion: ',NADD write(6,'(a)') '(This is effective only for spherical particles)' 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 C component particles is defined as f=(r_i+r_j)/d_{ij}, where (r_i,r_j) C are the radii of the two spheres circumscribing the particles i and j, C respectively, and d_{ij} is the separation distance of the two C particle-centers. When f', nL write(6,*) + 'Change nLp in gmm01f.par, recompile, then try again' stop endif if(nL.eq.1) idMie=1 C----------------------------------------------------------------------- C input the configuration and particle parameters for an snsemble C as mentioned above, one line for a particle and each line contains C 13 numbers: C idshp(i) - to specify particle shape C 0 for a sphere, -1 for a spheroid, -2 for a finite cylinder, ... C shp(1,i) - aspect ratio for a nonspherical particle C shp(2,i) and shp(3,i) - the azimuth and zenith angles of the C axis of rotational symmetry of a nonspherical particle C the rest of the line is the same as in the input data for the C codes gmm02f.f, gmm02s.f, or gmm02TrA.f, which includes 9 numbers: C x-, y-, z-coordinates of a particle-center, the radius of C an (equivalent) sphere in the same unit of the incident C wavelength, the real and imaginary parts of the refractive C index, the last three numbers are for a core-mantle sphere, C which are the radius, the real and imaginary parts of C refractive index of the core, for all other shapes of particles, C simply set these three to 0. C----------------------------------------------------------------------- do 1 i=1,nL read(2,*,err=10) idshp(i),(shp(j,i),j=1,3), + (r0(j,i),j=1,9) if(idshp(i).eq.0) then do j=1,3 shp(j,i)=0.d0 enddo endif temp=dabs(shp(1,i)-1.d0) temp0=1.d-7 if(idshp(i).eq.-1.and.temp.lt.temp0) then idshp(i)=0 do j=1,3 shp(j,i)=0.d0 enddo endif 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 if(irat.eq.1) then write(6,'(a,f7.3)') ' volume-equiv. xv: ', xv else write(6,'(a,f7.3)') ' surface-equiv. xs: ', xs endif write(6,'(/)') fileout='gmm03fq.out' fileout1=fileout if(idMie.eq.1) then write(6,*) + '*** Calculating only coherent scattering ***' write(6,*) + '*** No interaction included ****************' endif do i=1,nL x(i)=k*r0(4,i) if(r0(7,i).lt.1.d-12) then xc(i)=0.d0 refc(i)=dcmplx(0.d0,0.d0) else xc(i)=k*r0(7,i) endif enddo C ---------------------------------------------------------------------- C Calculating individual T-matrix for each component particles C For spherical particles, it only needs to calculate the "Mie" C scattering coefficients C (1) the subroutine "scoatabd.f" used here for calculating scattering C coefficients of homogeneous and core-mantle spheres is originally C the code "SCSERD.FOR" written by R.T. Wang and W.X. Wang C [see R.T. Wang and W.X. Wang, "Scattering by Arbitrarily Large C Homogeneous/Concentric Speres - Exact Theory with Use of New C Efficient Algorithms," in Proc. of the 1985 CRDC Scientific C Conference on Obscuration and Aerosol Research, R. Kohl, ed. C (Army CRDEC-SP-86019, Aberdeen, MD 1986), pp. 381-409], which C uses ratio method of Wang and van de Hulst in the calculation of C Riccati-Bessel functions [see Wang and van de Hulst, Appl. Opt. C 30, 106 (1991), Xu, Gustafson, Giovane, Blum, and Tehranian, C Phys. Rev. E 60, 2347 (1999)] C (2) for nonspherical particles, the individual-particle T-matrices C in the particle reference system are calculated by the subroutine C "tm0q.f", which is a part of Mishchenko's public domain code C "tmq.new.f" C ---------------------------------------------------------------------- nmax0=1 do i=1,nL do j1=1,nmp do j2=1,nmp tbar(i,1,1,j1,j2)=0.d0 tbar(i,1,2,j1,j2)=0.d0 tbar(i,2,1,j1,j2)=0.d0 tbar(i,2,2,j1,j2)=0.d0 enddo enddo if(i.eq.1) goto 12 do 121 j=i-1,1,-1 if(idshp(i).eq.idshp(j).and.shp(1,i).eq.shp(1,j)) then if(xc(i).eq.xc(j).and.refc(i).eq.refc(j)) then if(x(i).eq.x(j).and.ref(i).eq.ref(j)) then nmax(i)=nmax(j) uvmax(i)=uvmax(j) do n=1,nmax(j) do v=1,nmax(j) in0=n*(n+1) iv0=v*(v+1) tbar(i,1,1,in0,iv0)=tbar(j,1,1,in0,iv0) tbar(i,1,2,in0,iv0)=tbar(j,1,2,in0,iv0) tbar(i,2,1,in0,iv0)=tbar(j,2,1,in0,iv0) tbar(i,2,2,in0,iv0)=tbar(j,2,2,in0,iv0) do m=1,min(n,v) imn=m+in0 iuv=m+iv0 A=tbar(j,1,1,imn,iuv) tbar(i,1,1,imn,iuv)=A A=tbar(j,1,2,imn,iuv) tbar(i,1,2,imn,iuv)=A A=tbar(j,2,1,imn,iuv) tbar(i,2,1,imn,iuv)=A A=tbar(j,2,2,imn,iuv) tbar(i,2,2,imn,iuv)=A imn=-m+in0 iuv=-m+iv0 A=tbar(j,1,1,imn,iuv) tbar(i,1,1,imn,iuv)=A A=tbar(j,1,2,imn,iuv) tbar(i,1,2,imn,iuv)=A A=tbar(j,2,1,imn,iuv) tbar(i,2,1,imn,iuv)=A A=tbar(j,2,2,imn,iuv) tbar(i,2,2,imn,iuv)=A enddo enddo enddo c write(6,*) c + 'same type of particles: ',i,' and ',j goto 15 endif endif endif 121 continue 12 if(idshp(i).eq.0) then 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) ratio=0.d0 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,' write(6,*) ' recompile, then try again' stop endif uvmax(i)=nmax(i)*(nmax(i)+2) write(6,'(a,1x,i4)') + 'Actual single-particle expansion truncation:', + nmax(i) do j=1,uvmax(i) v=dsqrt(dble(j)) tbar(i,1,1,j,j)=an(v) tbar(i,2,2,j,j)=bn(v) enddo write(6,'(a,i3,a,f10.4)') 'sphere #',i, + ' individual size parameter: ',x(i) goto 15 endif RAT=irat DDELT=0.001D0 NDGS=2 do m=1,NPN6 do n=1,NPN4 do v=1,NPN4 RT11(m,n,v)=0.d0 RT12(m,n,v)=0.d0 RT21(m,n,v)=0.d0 RT22(m,n,v)=0.d0 IT11(m,n,v)=0.d0 IT12(m,n,v)=0.d0 IT21(m,n,v)=0.d0 IT22(m,n,v)=0.d0 enddo enddo enddo call tm0q(w,idshp(i),shp(1,i),r0(4,i),RAT, + r0(5,i),r0(6,i),DDELT,NDGS,nmax(i)) write(6,'(a,i3,a,f10.4)') 'particle #',i, + ' individual size parameter: ',x(i) write(6,*) 'individual T-matrix: ', i, ' ', 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,' write(6,*) ' recompile,then try again' stop endif if(nmax(i).gt.NPN1) then write(6,*) ' Parameter NPN1 too small, must be >', + nmax(i) write(6,*) ' Please change NPN1 in ampld.par.f,' write(6,*) ' recompile,then try again' stop endif uvmax(i)=nmax(i)*(nmax(i)+2) do n=1,nmax(i) in0=n*(n+1) do v=1,nmax(i) iv0=v*(v+1) nsmall=min(n,v) A=-ci**(v-n) imn=in0 iuv=iv0 B=dcmplx(RT11(1,n,v),IT11(1,n,v)) tbar(i,1,1,imn,iuv)=A*B B=dcmplx(RT12(1,n,v),IT12(1,n,v)) tbar(i,1,2,imn,iuv)=A*B B=dcmplx(RT21(1,n,v),IT21(1,n,v)) tbar(i,2,1,imn,iuv)=A*B B=dcmplx(RT22(1,n,v),IT22(1,n,v)) tbar(i,2,2,imn,iuv)=A*B do m=1,nsmall imn=m+in0 iuv=m+iv0 imn1=-m+in0 iuv1=-m+iv0 B=dcmplx(RT11(m+1,n,v),IT11(m+1,n,v)) tbar(i,1,1,imn,iuv)=A*B tbar(i,1,1,imn1,iuv1)=tbar(i,1,1,imn,iuv) B=dcmplx(RT12(m+1,n,v),IT12(m+1,n,v)) tbar(i,1,2,imn,iuv)=A*B tbar(i,1,2,imn1,iuv1)=-tbar(i,1,2,imn,iuv) B=dcmplx(RT21(m+1,n,v),IT21(m+1,n,v)) tbar(i,2,1,imn,iuv)=A*B tbar(i,2,1,imn1,iuv1)=-tbar(i,2,1,imn,iuv) B=dcmplx(RT22(m+1,n,v),IT22(m+1,n,v)) tbar(i,2,2,imn,iuv)=A*B tbar(i,2,2,imn1,iuv1)=tbar(i,2,2,imn,iuv) enddo enddo enddo 15 if(nmax(i).gt.nmax0) then nmax0=nmax(i) imax=i endif enddo write(6,*) 'maximum scattering order: ',imax,' ',nmax0 write(6,'(/)') write(6,*) 'input particle-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 calculating T-matrices of individual particles in their respective C specified orientations C-------------------------------------------------------------------- do 170 i=1,nL if(idshp(i).eq.0) goto 170 if(idalpha.eq.1) shp(2,i)=360.d0*ran1d(iseed) if(idbeta.eq.1) shp(3,i)=180.d0*ran1d(iseed) 170 enddo do 17 i=1,nL if(idshp(i).eq.0) goto 17 alph=shp(2,i)*pih/90.d0 beta=shp(3,i)*pih/90.d0 s=dabs(alph) t=dabs(beta) if(s.lt.1.d-10.and.t.lt.1.d-10) goto 17 if(i.eq.1) goto 172 do 171 j=i-1,1,-1 ii=idshp(i) ij=idshp(j) ca=shp(2,j) s=shp(1,i) t=shp(1,j) if(ii.eq.ij.and.s.eq.t) then if(x(i).eq.x(j).and.ref(i).eq.ref(j)) then sa=shp(2,i) ca=shp(2,j) sb=shp(3,i) cb=shp(3,j) if(sa.eq.ca.and.sb.eq.cb) then do in=1,uvmax(j) do iv=1,uvmax(j) tbar(i,1,1,in,iv)=tbar(j,1,1,in,iv) tbar(i,1,2,in,iv)=tbar(j,1,2,in,iv) tbar(i,2,1,in,iv)=tbar(j,2,1,in,iv) tbar(i,2,2,in,iv)=tbar(j,2,2,in,iv) enddo enddo goto 17 endif endif endif 171 continue 172 do in=1,uvmax(i) do iv=1,uvmax(i) tbar0(1,1,in,iv)=tbar(i,1,1,in,iv) tbar0(1,2,in,iv)=tbar(i,1,2,in,iv) tbar0(2,1,in,iv)=tbar(i,2,1,in,iv) tbar0(2,2,in,iv)=tbar(i,2,2,in,iv) enddo enddo sa=dsin(alph) ca=dcos(alph) cb=dcos(beta) A=dcmplx(ca,sa) n1=nmax(i) ekt(0)=1.d0 do m=1,2*n1 ekt(m)=A**m ekt(-m)=dconjg(ekt(m)) enddo call rotcoef(cb,n1) do m=-n1,n1 do u=-n1,n1 A=ekt(u-m) do 173 n=1,n1 if(iabs(m).gt.n) go to 173 in0=n*(n+1) imn=in0+m do 174 v=1,n1 if(iabs(u).gt.v) go to 173 iv0=v*(v+1) iuv=iv0+u nsmall=min(n,v) do ip=1,2 do iq=1,2 bar(ip,iq)=0.d0 enddo enddo do is=-nsmall,nsmall isn=in0+is isv=iv0+is t=dc(m,isn)*dc(u,isv) B=A*t do ip=1,2 do iq=1,2 Bj=tbar0(ip,iq,isn,isv) bar(ip,iq)=bar(ip,iq)+B*Bj enddo enddo enddo do ip=1,2 do iq=1,2 tbar(i,ip,iq,imn,iuv)=bar(ip,iq) enddo enddo 174 continue 173 continue enddo enddo 17 continue 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 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. C Phys. 139, 137 (1998) C call cofsrd(nmax0) call cofd0(nmax0) call cofnv0(nmax0) call gau0(nmax0) indpol=0 write(6,*) 'Starting Bi-CGSTAB solution process' write(6,*) 'Solving x-pol. inci. state' 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 and non-interacting C scattering coefficients for each component particles 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-2)=-A q0(i,imn-2)=A if(indpol.gt.1) then A=A*cin p0(i,imn-2)=A q0(i,imn-2)=-A endif if(idshp(i).eq.0) then B=tbar(i,1,1,imn,imn) as(i,imn)=B*A as(i,imn-2)=B*p0(i,imn-2) B=tbar(i,2,2,imn,imn) bs(i,imn)=B*A bs(i,imn-2)=B*q0(i,imn-2) else do iuv=1,uvmax(i) B=tbar(i,1,1,iuv,imn) B=B+tbar(i,1,2,iuv,imn) A0=B*A as(i,iuv)=as(i,iuv)+A0 B=tbar(i,1,1,iuv,imn-2) A0=B*p0(i,imn-2) B=tbar(i,1,2,iuv,imn-2) A0=A0+B*q0(i,imn-2) as(i,iuv)=as(i,iuv)+A0 B=tbar(i,2,1,iuv,imn) B=B+tbar(i,2,2,iuv,imn) A0=B*A bs(i,iuv)=bs(i,iuv)+A0 B=tbar(i,2,1,iuv,imn-2) A0=B*p0(i,imn-2) B=tbar(i,2,2,iuv,imn-2) A0=A0+B*q0(i,imn-2) bs(i,iuv)=bs(i,iuv)+A0 enddo endif 21 continue if(idshp(i).eq.0) then do n=1,nmax(i) imn=n*n+n+1 p0(i,imn)=as(i,imn) q0(i,imn)=bs(i,imn) p0(i,imn-2)=as(i,imn-2) q0(i,imn-2)=bs(i,imn-2) enddo else do iuv=1,uvmax(i) p0(i,iuv)=as(i,iuv) q0(i,iuv)=bs(i,iuv) enddo endif 20 continue C----------------------------------------------------------------- c begins iteration process to solve the interaction equations for c partial interacting scatteing coefficients using BI-CGSTAB [see c H.A. van der Vorst, SIAM, J. Sci. Stat. Comput. 13, 631 (1992); c Gutknecht, SIAM J. Sci. Comput. 14, 1020 (1993)] C--------------------------------------------------------------- if(idMie.eq.1) goto 490 if(nL.eq.1) goto 490 do i=1,nL ind(i)=0 c0i(i)=0.d0 if(idshp(i).eq.0) then do n=1,nmax(i) j=n*n+n+1 c0i(i)=c0i(i)+p0(i,j)*dconjg(p0(i,j)) c0i(i)=c0i(i)+q0(i,j)*dconjg(q0(i,j)) c0i(i)=c0i(i)+p0(i,j-2)*dconjg(p0(i,j-2)) c0i(i)=c0i(i)+q0(i,j-2)*dconjg(q0(i,j-2)) enddo else do j=1,uvmax(i) c0i(i)=c0i(i)+p0(i,j)*dconjg(p0(i,j)) c0i(i)=c0i(i)+q0(i,j)*dconjg(q0(i,j)) enddo endif enddo icase=1 niter=1 call transs(k,nL,r0,nmax0,nmax,uvmax,fint,ek, + besj,besy,drot,as,bs,as00,bs00,ind,icase) do i=1,nL c1i(i)=0.d0 if(idshp(i).eq.0) then do imn=1,uvmax(i) s2x=tbar(i,1,1,imn,imn) as1(i,imn)=-s2x*as00(i,imn) s4x=tbar(i,2,2,imn,imn) bs1(i,imn)=-s4x*bs00(i,imn) c1i(i)=c1i(i)+as1(i,imn)*dconjg(as1(i,imn)) c1i(i)=c1i(i)+bs1(i,imn)*dconjg(bs1(i,imn)) enddo else if(shp(2,i).eq.0.d0.and.shp(3,i).eq.0.d0) then do imn=1,uvmax(i) n=dsqrt(dble(imn)) m=imn-n*(n+1) s3y=0.d0 s1y=0.d0 do v=1,nmax(i) if(iabs(m).le.v) then iuv=m+v*(v+1) s2x=tbar(i,1,1,imn,iuv) s3y=s3y+s2x*as00(i,iuv) s2x=tbar(i,1,2,imn,iuv) s3y=s3y+s2x*bs00(i,iuv) s4x=tbar(i,2,1,imn,iuv) s1y=s1y+s4x*as00(i,iuv) s4x=tbar(i,2,2,imn,iuv) s1y=s1y+s4x*bs00(i,iuv) endif enddo as1(i,imn)=-s3y bs1(i,imn)=-s1y c1i(i)=c1i(i)+s3y*dconjg(s3y) c1i(i)=c1i(i)+s1y*dconjg(s1y) enddo else do imn=1,uvmax(i) s3y=0.d0 s1y=0.d0 do iuv=1,uvmax(i) s2x=tbar(i,1,1,imn,iuv) s3y=s3y+s2x*as00(i,iuv) s2x=tbar(i,1,2,imn,iuv) s3y=s3y+s2x*bs00(i,iuv) s4x=tbar(i,2,1,imn,iuv) s1y=s1y+s4x*as00(i,iuv) s4x=tbar(i,2,2,imn,iuv) s1y=s1y+s4x*bs00(i,iuv) enddo as1(i,imn)=-s3y bs1(i,imn)=-s1y c1i(i)=c1i(i)+s3y*dconjg(s3y) c1i(i)=c1i(i)+s1y*dconjg(s1y) enddo endif endif enddo temp=0.d0 do i=1,nL cext0=c1i(i)/c0i(i) if(cext0.lt.small) ind(i)=1 if(cext0.gt.temp) temp=cext0 enddo if(temp.lt.small) goto 490 A0=0.d0 do i=1,nL if(ind(i).gt.0) goto 611 do imn=1,uvmax(i) asp(i,imn)=as1(i,imn) bsp(i,imn)=bs1(i,imn) as0(i,imn)=as1(i,imn) bs0(i,imn)=bs1(i,imn) A0=A0+as1(i,imn)*as1(i,imn) A0=A0+bs1(i,imn)*bs1(i,imn) enddo 611 continue enddo 62 call transs(k,nL,r0,nmax0,nmax,uvmax,fint,ek, + besj,besy,drot,asp,bsp,as00,bs00,ind,icase) do 621 i=1,nL if(ind(i).gt.0) goto 621 if(idshp(i).eq.0) then do imn=1,uvmax(i) s2x=tbar(i,1,1,imn,imn) asv(i,imn)=s2x*as00(i,imn) s4x=tbar(i,2,2,imn,imn) bsv(i,imn)=s4x*bs00(i,imn) asv(i,imn)=asv(i,imn)+asp(i,imn) bsv(i,imn)=bsv(i,imn)+bsp(i,imn) enddo else if(shp(2,i).eq.0.d0.and.shp(3,i).eq.0.d0) then do imn=1,uvmax(i) n=dsqrt(dble(imn)) m=imn-n*(n+1) s3y=0.d0 s1y=0.d0 do v=1,nmax(i) if(iabs(m).le.v) then iuv=m+v*(v+1) s2x=tbar(i,1,1,imn,iuv) s3y=s3y+s2x*as00(i,iuv) s2x=tbar(i,1,2,imn,iuv) s3y=s3y+s2x*bs00(i,iuv) s4x=tbar(i,2,1,imn,iuv) s1y=s1y+s4x*as00(i,iuv) s4x=tbar(i,2,2,imn,iuv) s1y=s1y+s4x*bs00(i,iuv) endif enddo asv(i,imn)=s3y bsv(i,imn)=s1y asv(i,imn)=asv(i,imn)+asp(i,imn) bsv(i,imn)=bsv(i,imn)+bsp(i,imn) enddo else do imn=1,uvmax(i) s3y=0.d0 s1y=0.d0 do iuv=1,uvmax(i) s2x=tbar(i,1,1,imn,iuv) s3y=s3y+s2x*as00(i,iuv) s2x=tbar(i,1,2,imn,iuv) s3y=s3y+s2x*bs00(i,iuv) s4x=tbar(i,2,1,imn,iuv) s1y=s1y+s4x*as00(i,iuv) s4x=tbar(i,2,2,imn,iuv) s1y=s1y+s4x*bs00(i,iuv) enddo asv(i,imn)=s3y bsv(i,imn)=s1y asv(i,imn)=asv(i,imn)+asp(i,imn) bsv(i,imn)=bsv(i,imn)+bsp(i,imn) enddo endif endif 621 continue A=0.d0 do 622 i=1,nL if(ind(i).gt.0) goto 622 do imn=1,uvmax(i) A=A+asv(i,imn)*as1(i,imn) A=A+bsv(i,imn)*bs1(i,imn) enddo 622 continue Aj=A0/A do 623 i=1,nL if(ind(i).gt.0) goto 623 do imn=1,uvmax(i) asc(i,imn)=as0(i,imn)-Aj*asv(i,imn) bsc(i,imn)=bs0(i,imn)-Aj*bsv(i,imn) enddo 623 continue call transs(k,nL,r0,nmax0,nmax,uvmax,fint,ek, + besj,besy,drot,asc,bsc,as00,bs00,ind,icase) do 624 i=1,nL if(ind(i).gt.0) goto 624 if(idshp(i).eq.0) then do imn=1,uvmax(i) s2x=tbar(i,1,1,imn,imn) ast(i,imn)=s2x*as00(i,imn) s4x=tbar(i,2,2,imn,imn) bst(i,imn)=s4x*bs00(i,imn) ast(i,imn)=ast(i,imn)+asc(i,imn) bst(i,imn)=bst(i,imn)+bsc(i,imn) enddo else if(shp(2,i).eq.0.d0.and.shp(3,i).eq.0.d0) then do imn=1,uvmax(i) n=dsqrt(dble(imn)) m=imn-n*(n+1) s3y=0.d0 s1y=0.d0 do v=1,nmax(i) if(iabs(m).le.v) then iuv=m+v*(v+1) s2x=tbar(i,1,1,imn,iuv) s3y=s3y+s2x*as00(i,iuv) s2x=tbar(i,1,2,imn,iuv) s3y=s3y+s2x*bs00(i,iuv) s4x=tbar(i,2,1,imn,iuv) s1y=s1y+s4x*as00(i,iuv) s4x=tbar(i,2,2,imn,iuv) s1y=s1y+s4x*bs00(i,iuv) endif enddo ast(i,imn)=s3y bst(i,imn)=s1y ast(i,imn)=ast(i,imn)+asc(i,imn) bst(i,imn)=bst(i,imn)+bsc(i,imn) enddo else do imn=1,uvmax(i) s3y=0.d0 s1y=0.d0 do iuv=1,uvmax(i) s2x=tbar(i,1,1,imn,iuv) s3y=s3y+s2x*as00(i,iuv) s2x=tbar(i,1,2,imn,iuv) s3y=s3y+s2x*bs00(i,iuv) s4x=tbar(i,2,1,imn,iuv) s1y=s1y+s4x*as00(i,iuv) s4x=tbar(i,2,2,imn,iuv) s1y=s1y+s4x*bs00(i,iuv) enddo ast(i,imn)=s3y bst(i,imn)=s1y ast(i,imn)=ast(i,imn)+asc(i,imn) bst(i,imn)=bst(i,imn)+bsc(i,imn) enddo endif endif 624 continue A2=0.d0 B2=0.d0 do 625 i=1,nL if(ind(i).gt.0) goto 625 do imn=1,uvmax(i) A2=A2+ast(i,imn)*asc(i,imn) A2=A2+bst(i,imn)*bsc(i,imn) B2=B2+ast(i,imn)*ast(i,imn) B2=B2+bst(i,imn)*bst(i,imn) enddo 625 continue Bj=A2/B2 do 626 i=1,nL if(ind(i).gt.0) goto 626 c1i(i)=0.d0 do imn=1,uvmax(i) Aj2=Aj*asp(i,imn)+Bj*asc(i,imn) Bj2=Aj*bsp(i,imn)+Bj*bsc(i,imn) c1i(i)=c1i(i)+Aj2*dconjg(Aj2) c1i(i)=c1i(i)+Bj2*dconjg(Bj2) as(i,imn)=as(i,imn)+Aj2 bs(i,imn)=bs(i,imn)+Bj2 enddo 626 continue cext0=0.d0 cext1=0.d0 do 627 i=1,nL if(ind(i).gt.0) goto 627 cext0=cext0+c0i(i) cext1=cext1+c1i(i) 627 continue temp=cext1/cext0 if(temp.lt.small) goto 490 if(niter.gt.MXINT) then write(6,*) 'Caution:' write(6,*) '****** Maximum iterations exceeded ******' write(6,*) '****** Results may be inaccurate!! ******' goto 490 endif if(niter.eq.1.or.niter/50*50.eq.niter) + write(6,'(a,i5,2x,e15.7)') 'iteration # ',niter,temp do 628 i=1,nL if(ind(i).gt.0) goto 628 do imn=1,uvmax(i) as0(i,imn)=asc(i,imn)-Bj*ast(i,imn) bs0(i,imn)=bsc(i,imn)-Bj*bst(i,imn) enddo 628 continue A2=0.d0 do 629 i=1,nL ast(i,1)=0.d0 if(ind(i).gt.0) goto 629 do imn=1,uvmax(i) ast(i,1)=ast(i,1)+as0(i,imn)*as1(i,imn) ast(i,1)=ast(i,1)+bs0(i,imn)*bs1(i,imn) enddo A2=A2+ast(i,1) 629 continue B0=A2/A0 B0=B0*Aj/Bj do 630 i=1,nL if(ind(i).gt.0) goto 630 do imn=1,uvmax(i) asp(i,imn)=as0(i,imn)+B0*(asp(i,imn)-Bj*asv(i,imn)) bsp(i,imn)=bs0(i,imn)+B0*(bsp(i,imn)-Bj*bsv(i,imn)) enddo 630 continue A0=0.d0 do 631 i=1,nL if(ind(i).gt.0) goto 631 cext0=c1i(i)/c0i(i) if(cext0.lt.small) then ind(i)=1 goto 631 endif A0=A0+ast(i,1) 631 continue niter=niter+1 goto 62 490 indpol=indpol+2 if(indpol.lt.3) then do i=1,nL do j=1,uvmax(i) asx(i,j)=as(i,j) bsx(i,j)=bs(i,j) enddo enddo write(6,*) 'Solving y-pol. inci. state' goto 18 endif C---------------------------------------------------------------- c computing total and differential scattering cross sections and c efficiencies for radiation pressure, using the formulas given c by Xu, Physics Letters A 249, 30 (1998) C---------------------------------------------------------------- do i=1,nL ind(i)=0 enddo icase=2 do 700 indpol=0,2,2 if(indpol.lt.1) then call transs(k,nL,r0,nmax0,nmax,uvmax,fint,ek, + besj,besy,drot,asx,bsx,as1,bs1,ind,icase) else call transs(k,nL,r0,nmax0,nmax,uvmax,fint,ek, + besj,besy,drot,as,bs,as1,bs1,ind,icase) endif do i=1,nL do imn=1,uvmax(i) if(indpol.lt.1) then at(imn)=asx(i,imn)+as1(i,imn) bt(imn)=bsx(i,imn)+bs1(i,imn) else at(imn)=as(i,imn)+as1(i,imn) bt(imn)=bs(i,imn)+bs1(i,imn) endif enddo do n=1,nmax(i) n1=n+1 n2=2*n rn=1.0d0/dble(n*n1) p=fnr(n)*fnr(n+2)/fnr(n2+1)/fnr(n2+3)/dble(n1) t=fnr(n-1)*fnr(n+1)/fnr(n2-1)/fnr(n2+1)/dble(n) sc=0.d0 temp=0.d0 do m=-n,n iL=n*(n+1)+m if(indpol.lt.1) then sc=sc+dble(dconjg(asx(i,iL))*at(iL)) sc=sc+dble(dconjg(bsx(i,iL))*bt(iL)) else sc=sc+dble(dconjg(as(i,iL))*at(iL)) sc=sc+dble(dconjg(bs(i,iL))*bt(iL)) endif rm=dble(m)*rn A0=rm*bt(iL) B0=rm*at(iL) if(n.eq.nmax(i)) goto 51 u=(n+1)*(n+2)+m fnp=fnr(n+m+1)*fnr(n-m+1)*p A0=A0+fnp*at(u) B0=B0+fnp*bt(u) 51 if(n.eq.1.or.iabs(m).gt.n-1) goto 52 u=(n-1)*n+m fn=fnr(n+m)*fnr(n-m)*t A0=A0+fn*at(u) B0=B0+fn*bt(u) 52 if(indpol.lt.1) then temp=temp+dble(dconjg(asx(i,iL))*A0) temp=temp+dble(dconjg(bsx(i,iL))*B0) else temp=temp+dble(dconjg(as(i,iL))*A0) temp=temp+dble(dconjg(bs(i,iL))*B0) endif enddo if(indpol.lt.1) then cscaxi(i)=cscaxi(i)+sc cscax=cscax+sc cprxi(i)=cprxi(i)+temp cprx=cprx+temp else cscayi(i)=cscayi(i)+sc cscay=cscay+sc cpryi(i)=cpryi(i)+temp cpry=cpry+temp endif enddo enddo 700 continue C------------------------------------------------------------------- c computing total and differential extinction cross-sections c [see, for example, Xu, Appl. Opt. 36, 9496 (1997)] C------------------------------------------------------------------- do j=1,nL cz=dcos(k*r0(3,j)) sz=dsin(k*r0(3,j)) cmz=dcmplx(cz,-sz) A=0.d0 B=0.d0 Aj=0.d0 Bj=0.d0 do n=1,nmax(j) rn=fnr(2*n+1) m0=n*n+n+1 u0=n*n+n-1 A=A+rn*(asx(j,m0)+bsx(j,m0)) B=B+rn*(asx(j,u0)-bsx(j,u0)) Aj=Aj+rn*(as(j,m0)+bs(j,m0)) Bj=Bj+rn*(as(j,u0)-bs(j,u0)) enddo cextxi(j)=cextxi(j)+dble((A-B)*cmz) cextx=cextx+dble((A-B)*cmz) cextyi(j)=cextyi(j)-dimag((Aj+Bj)*cmz) cexty=cexty-dimag((Aj+Bj)*cmz) enddo C---------------------------------------------------------------------- c computing two-dimensional amplitude scattering matrix and c scattering Mueller matrix of an ensemble as a whole in an array of c (nang+1)*(npng+1) for both front and rear hemispheres c (nang and npng are determined by the input parameters sang and c pang, respectively) c Note that the definition used here for the amplitude scattering c matrix is slightly different from that given by van de Hulst and by c Bohren & Huffman. The incident field components refer to the x-z c plane instead of the scattering plane defined by the z axis and c the scattering direction. C---------------------------------------------------------------------- write(6,'(/)') write(6,*) ' ********************************************' write(6,*) ' Output files for the amplitude scattering' write(6,*) ' matrix are amp.forw.out and amp.back.out' write(6,*) ' and for the scattering Mueller matrix are' write(6,*) ' mue.forw.out and mue.back.out' write(6,*) ' ********************************************' write(6,'(/)') open(19,file='amp.forw.out',status='unknown') write(19,'(a)') + 'amp.forw.out (amplitude scattering matrix for the' write(19,'(a)') + ' hemisphere of forward scattering directions)' write(19,'(a,f7.4,3x,a,a20)') + 'wavelength: ',w,'input filename: ',FLNAME write(19,'(a)') + 'particle#,x,y,z,r,complex refractive index:' do i=1,nL,nL write(19,'(i5,6f11.4)') i,r0(1,i),r0(2,i), + r0(3,i),r0(4,i),r0(5,i),r0(6,i) enddo write(19,'(a)') 'for each scattering angle theta ', + 'in the range of [0,90] (degrees)' write(19,'(a)') + 'azimuth angle, s2x(complex), s3y(complex)' write(19,'(a)') + ' s4x(complex), s1y(complex)' open(20,file='amp.back.out',status='unknown') write(20,'(a)') + 'amp.forw.out (amplitude scattering matrix for the' write(20,'(a)') + ' hemisphere of backward scattering directions)' write(20,'(a,f7.4,3x,a,a20)') + 'wavelength: ',w,'input filename: ',FLNAME write(20,'(a)') + 'particle#,x,y,z,r,complex refractive index:' do i=1,nL,nL write(20,'(i5,6f11.4)') i,r0(1,i),r0(2,i), + r0(3,i),r0(4,i),r0(5,i),r0(6,i) enddo write(20,'(a)') 'for each scattering angle theta ', + 'in the range of [90,180] (degrees)' write(20,'(a)') + 'azimuth angle (phi), s2x(complex), s3y(complex)' write(20,'(a)') + ' s4x(complex), s1y(complex)' open(21,file='mue.forw.out',status='unknown') write(21,'(a)') 'mue.forw.out' write(21,'(a)') 'Mueller matrix for the hemisphere of forward' write(21,'(a)') 'scattering directions and for a single fixed' write(21,'(a)') 'orientation' write(21,'(a16,a20)') 'Input filename: ',FLNAME write(21,'(a30,f9.3)') + 'interval of scattering angle: ',sang write(21,'(a30,f9.3)') + 'interval of azimuth angle: ',pang open(22,file='mue.back.out',status='unknown') write(22,'(a)') 'mue.back.out' write(22,'(a)') 'Mueller matrix for the hemisphere of backward' write(22,'(a)') 'scattering directions and for a single fixed' write(22,'(a)') 'orientation' write(22,'(a16,a20)') 'Input filename: ',FLNAME write(22,'(a30,f9.3)') + 'interval of scattering angle: ',sang write(22,'(a30,f9.3)') + 'interval of azimuth angle: ',pang do i=1,nang iang=2*nang-i dang(i)=sang*dble(i-1) dang(iang)=180.0d0-dang(i) write(19,'(a18,f8.3)') 'theta in degrees: ',dang(i) write(20,'(a18,f8.3)') 'theta in degrees: ',dang(iang) write(21,'(a18,f8.3)') 'theta in degrees: ',dang(i) write(22,'(a18,f8.3)') 'theta in degrees: ',dang(iang) theta=dang(i)*pione/180.0d0 xt=dcos(theta) st=dsin(theta) call tipitaud(nmax0,xt) do jc=1,npng+1 phijc=pang*dble(jc-1) azphi=pih*phijc/90.0d0 sphi=dsin(azphi) cphi=dcos(azphi) do imn=1,nmp at(imn)=0.d0 bt(imn)=0.d0 atj(imn)=0.d0 btj(imn)=0.d0 aty(imn)=0.d0 bty(imn)=0.d0 atyj(imn)=0.d0 btyj(imn)=0.d0 enddo do 31 j=1,nL sb=r0(1,j)*cphi+r0(2,j)*sphi sb=sb*st cb=r0(3,j)*xt cz=k*(sb+cb) sz=k*(sb-cb) A=dcmplx(dcos(cz),-dsin(cz)) B=dcmplx(dcos(sz),-dsin(sz)) do 32 n=1,nmax(j) n1=n*(n+1) do 33 m=-n,n imn=m+n1 if(idMie.gt.0) then if(idshp(j).eq.0) then if(iabs(m).ne.1) goto 33 endif endif at(imn)=at(imn)+A*asx(j,imn) bt(imn)=bt(imn)+A*bsx(j,imn) atj(imn)=atj(imn)+B*asx(j,imn) btj(imn)=btj(imn)+B*bsx(j,imn) aty(imn)=aty(imn)+A*as(j,imn) bty(imn)=bty(imn)+A*bs(j,imn) atyj(imn)=atyj(imn)+B*as(j,imn) btyj(imn)=btyj(imn)+B*bs(j,imn) 33 continue 32 continue 31 continue s2x=0.d0 s4x=0.d0 s2xj=0.d0 s4xj=0.d0 s3y=0.d0 s1y=0.d0 s3yj=0.d0 s1yj=0.d0 A=0.d0 B=0.d0 Aj=0.d0 Bj=0.d0 A2=0.d0 B2=0.d0 Aj2=0.d0 Bj2=0.d0 do j=1,nmax0 imn=(j-1)*(j+2)/2+1 u=j*j+j A=A+at(u)*tau(imn) B=B+bt(u)*tau(imn) A2=A2+aty(u)*tau(imn) B2=B2+bty(u)*tau(imn) t=(-1)**(j+1) Aj=Aj+atj(u)*tau(imn)*t Bj=Bj+btj(u)*tau(imn)*t Aj2=Aj2+atyj(u)*tau(imn)*t Bj2=Bj2+btyj(u)*tau(imn)*t enddo s2x=s2x+A s4x=s4x+B*dcmplx(0.d0,-1.d0) s3y=s3y-A2 s1y=s1y+B2*dcmplx(0.d0,1.d0) s2xj=s2xj+Aj s4xj=s4xj+Bj*dcmplx(0.d0,-1.d0) s3yj=s3yj-Aj2 s1yj=s1yj+Bj2*dcmplx(0.d0,1.d0) rm=1.d0 do 302 m=1,nmax0 A=0.d0 B=0.d0 A2=0.d0 B2=0.d0 Aj=0.d0 Bj=0.d0 Aj2=0.d0 Bj2=0.d0 do jm=1,2 Ay(jm)=0.d0 By(jm)=0.d0 Ayj(jm)=0.d0 Byj(jm)=0.d0 enddo rm=-rm do 303 j=m,nmax0 imn=(j-1)*(j+2)/2+m+1 u=j*j+j+m v=u-2*m A0=at(u)*tau(imn)+bt(u)*pi(imn) B0=rm*(at(v)*tau(imn)-bt(v)*pi(imn)) A=A+A0+B0 A2=A2+A0-B0 A0=aty(u)*tau(imn)+bty(u)*pi(imn) B0=rm*(aty(v)*tau(imn)-bty(v)*pi(imn)) Ay(1)=Ay(1)+A0+B0 Ay(2)=Ay(2)+A0-B0 A0=at(u)*pi(imn)+bt(u)*tau(imn) B0=rm*(at(v)*pi(imn)-bt(v)*tau(imn)) B=B+A0-B0 B2=B2+A0+B0 A0=aty(u)*pi(imn)+bty(u)*tau(imn) B0=rm*(aty(v)*pi(imn)-bty(v)*tau(imn)) By(1)=By(1)+A0-B0 By(2)=By(2)+A0+B0 t=(-1)**(j+m+1) p=-t A0=atj(u)*tau(imn)*t+btj(u)*pi(imn)*p B0=rm*(atj(v)*tau(imn)*t-btj(v)*pi(imn)*p) Aj=Aj+A0+B0 Aj2=Aj2+A0-B0 A0=atyj(u)*tau(imn)*t+btyj(u)*pi(imn)*p B0=rm*(atyj(v)*tau(imn)*t-btyj(v)*pi(imn)*p) Ayj(1)=Ayj(1)+A0+B0 Ayj(2)=Ayj(2)+A0-B0 A0=atj(u)*pi(imn)*p+btj(u)*tau(imn)*t B0=rm*(atj(v)*pi(imn)*p-btj(v)*tau(imn)*t) Bj=Bj+A0-B0 Bj2=Bj2+A0+B0 A0=atyj(u)*pi(imn)*p+btyj(u)*tau(imn)*t B0=rm*(atyj(v)*pi(imn)*p-btyj(v)*tau(imn)*t) Byj(1)=Byj(1)+A0-B0 Byj(2)=Byj(2)+A0+B0 303 continue temp=dble(m)*azphi sb=dsin(temp) cb=dcos(temp) s2x=s2x+A*cb+A2*dcmplx(0.d0,1.d0)*sb s4x=s4x+B*dcmplx(0.d0,-1.d0)*cb+B2*sb s3y=s3y-Ay(1)*cb+Ay(2)*dcmplx(0.d0,-1.d0)*sb s1y=s1y+By(1)*dcmplx(0.d0,1.d0)*cb-By(2)*sb s2xj=s2xj+Aj*cb+Aj2*dcmplx(0.d0,1.d0)*sb s4xj=s4xj+Bj*dcmplx(0.d0,-1.d0)*cb+Bj2*sb s3yj=s3yj-Ayj(1)*cb+Ayj(2)*dcmplx(0.d0,-1.d0)*sb s1yj=s1yj+Byj(1)*dcmplx(0.d0,1.d0)*cb-Byj(2)*sb 302 continue write(19,'(f8.2,4e15.6)') phijc, + dble(s2x),dimag(s2x),dble(s3y),dimag(s3y) write(19,'(8x,4e15.6)') + dble(s4x),dimag(s4x),dble(s1y),dimag(s1y) write(20,'(f8.2,4e15.6)') phijc, + dble(s2xj),dimag(s2xj),dble(s3yj),dimag(s3yj) write(20,'(8x,4e15.6)') + dble(s4xj),dimag(s4xj),dble(s1yj),dimag(s1yj) if(jc.eq.1) then i22(i)=cdabs(s2x)*cdabs(s2x) i21(i)=cdabs(s4x)*cdabs(s4x) i11(i)=cdabs(s1y)*cdabs(s1y) i12(i)=cdabs(s3y)*cdabs(s3y) i22(iang)=cdabs(s2xj)*cdabs(s2xj) i21(iang)=cdabs(s4xj)*cdabs(s4xj) i11(iang)=cdabs(s1yj)*cdabs(s1yj) i12(iang)=cdabs(s3yj)*cdabs(s3yj) if(i.eq.1) then cbakx=cdabs(s2xj)*cdabs(s2xj) cbaky=cdabs(s1yj)*cdabs(s1yj) endif endif call mueller(s1y,s2x,s3y,s4x,smue) write(21,'(f7.1,4e16.7)') + phijc,smue(1,1),smue(1,2),smue(1,3),smue(1,4) write(21,'(7x,4e16.7)') + smue(2,1),smue(2,2),smue(2,3),smue(2,4) write(21,'(7x,4e16.7)') + smue(3,1),smue(3,2),smue(3,3),smue(3,4) write(21,'(7x,4e16.7)') + smue(4,1),smue(4,2),smue(4,3),smue(4,4) call mueller(s1yj,s2xj,s3yj,s4xj,smue) write(22,'(f7.1,4e16.7)') + phijc,smue(1,1),smue(1,2),smue(1,3),smue(1,4) write(22,'(7x,4e16.7)') + smue(2,1),smue(2,2),smue(2,3),smue(2,4) write(22,'(7x,4e16.7)') + smue(3,1),smue(3,2),smue(3,3),smue(3,4) write(22,'(7x,4e16.7)') + smue(4,1),smue(4,2),smue(4,3),smue(4,4) enddo enddo close(19) close(20) close(21) close(22) 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 cz=k*k 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=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=(cprx+cpry)/(cscax+cscay) assym0=0.5d0*(cprx/cscax+cpry/cscay) cbakx=2.0d0*twopi*cbakx/cz cbaky=2.0d0*twopi*cbaky/cz cbak=0.5d0*(cbakx+cbaky) c write(6,'(5x,2e15.6)') assym,assym0 c write(6,'(6x,a4,9x,a4,9x,a4,9x,a4,9x,a3,8x,a12)') c + 'Cext','Cabs','Csca','Cbak','Cpr','' c write(6,'(2x,6e13.5)') cext,cabs,csca,cbak,cext-cpr,assym cscax=0.d0 cscay=0.d0 cextx=0.d0 cexty=0.d0 cprx=0.d0 cpry=0.d0 do i=1,nL cscax=cscax+cscaxi(i) cscay=cscay+cscayi(i) cextx=cextx+cextxi(i) cexty=cexty+cextyi(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 c write(6,'(2x,6e13.5)') cext,cabs,csca,cbak,cext-cpr,assym do i=1,nL 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) c write(6,'(i5,5e15.6)') i,cexti(i),cabsi(i), c + cscai(i),cpri(i),assymi(i) enddo c write(6,'(/,a)') 'efficiencies for radiation pressure' c write(6,'(5x,6e13.5)') assym,assym0 c write(6,'(5x,6e13.5)') assym,cext-cpr,assymx,cextx-cprx, c + assymy,cexty-cpry c do i=1,nL c write(6,'(i5,6e13.5)') i,assymi(i),cpri(i), c + assymxi(i),cprxi(i),assymyi(i),cpryi(i) c enddo 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,*) '(a single fixed orientation)' 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) if(irat.eq.1) then cz=pione*gcvr*gcvr else cz=pione*gcsr*gcsr endif 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) 222 format(1x,a1,6e13.5) 221 format(6x,a5,8x,a5,8x,a5,8x,a5,8x,a4,5x,a12) if(irat.eq.1) then write(6,221) + 'Qextv','Qabsv','Qscav','Qbakv','Qprv','' else write(6,221) + 'Qexts','Qabss','Qscas','Qbaks','Qprs','' endif 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 open(12,file=fileout,status='unknown') write(12,*) 'Numerical solutions from the gmm03fq.f code' write(12,*) 'for a single or ensembles of variouly shaped' write(12,*) 'small particles in a single, fixed orientation' if(irat.eq.1) then write(12,'(a20,a16,a18,a4,f8.3)') fileout, + '--- input file: ',FLNAME,' xv:',xv else write(12,'(a20,a16,a18,a4,f8.3)') fileout, + '--- input file: ',FLNAME,' xs:',xs endif 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 if(irat.eq.1) then write(12,221) + 'Qextv','Qabsv','Qscav','Qbakv','Qprv','' else write(12,221) + 'Qexts','Qabss','Qscas','Qbaks','Qprs','' endif 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,'(/)') 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) STOP 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 + -fnr(n-k)*fnr(n+k+1)*sbe2*dkn1 + -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 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) 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 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 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) 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 C ------------------------------------------------------------------- C C the following is a part of the public domain code "tmq.new.f" C available online at URL: www.giss.nasa.gov/~crmin. For detailed C information on the formulation and algorithms, please see the C original code by Mishchenko C C ------------------------------------------------------------------- subroutine tm0q(LAM0,NP,EPS0,AXMAX,RAT,MRR0,MRI0, + DDELT,NDGS,NMAX) IMPLICIT double precision (A-H,O-Z) INCLUDE 'tmq.par.f' REAL*16 LAM,MRR,MRI,X(NPNG2),W(NPNG2),S(NPNG2),SS(NPNG2), * AN(NPN1),R(NPNG2),DR(NPNG2),PPI,PIR,PII,P,EPS,A, * DDR(NPNG2),DRR(NPNG2),DRI(NPNG2),ANN(NPN1,NPN1) double precision LAM0,MRR0,MRI0, & XG(1000),WG(1000),TR1(NPN2,NPN2),TI1(NPN2,NPN2), & ALPH1(NPL),ALPH2(NPL),ALPH3(NPL),ALPH4(NPL),BET1(NPL), & BET2(NPL),XG1(2000),WG1(2000), & AL1(NPL),AL2(NPL),AL3(NPL),AL4(NPL),BE1(NPL),BE2(NPL) double precision & RT11(NPN6,NPN4,NPN4),RT12(NPN6,NPN4,NPN4), & RT21(NPN6,NPN4,NPN4),RT22(NPN6,NPN4,NPN4), & IT11(NPN6,NPN4,NPN4),IT12(NPN6,NPN4,NPN4), & IT21(NPN6,NPN4,NPN4),IT22(NPN6,NPN4,NPN4) COMMON /CT/ TR1,TI1 COMMON /TMAT/ RT11,RT12,RT21,RT22,IT11,IT12,IT21,IT22 COMMON /CHOICE/ ICHOICE LAM=LAM0 MRR=MRR0 MRI=MRI0 EPS=EPS0 P=DACOS(-1.d0) PIN=P C OPEN FILES ******************************************************* c OPEN (6,FILE='test') c OPEN (10,FILE='tmatr.write') C INPUT DATA ******************************************************** NDISTR=4 NPNAX=1 B=1D-1 GAM=0.5D0 NKMAX=-1 NPNA=19 NDGS=2 NPNAX=1 ICHOICE=2 NCHECK=0 IF (NP.EQ.-1.OR.NP.EQ.-2) NCHECK=1 IF (NP.GT.0.AND.(-1)**NP.EQ.1) NCHECK=1 WRITE (6,5454) ICHOICE,NCHECK 5454 FORMAT ('ICHOICE=',I1,' NCHECK=',I1) DAX=AXMAX/NPNAX DLAM=LAM DMRR=MRR DMRI=MRI DEPS=EPS IF (DABS(RAT-1D0).GT.1D-8.AND.NP.EQ.-1) CALL SAREA (DEPS,RAT) IF (DABS(RAT-1D0).GT.1D-8.AND.NP.GE.0) CALL SURFCH(NP,DEPS,RAT) IF (DABS(RAT-1D0).GT.1D-8.AND.NP.EQ.-2) CALL SAREAC (DEPS,RAT) C PRINT 8000, RAT 8000 FORMAT ('RAT=',F8.6) IF(NP.EQ.-1.AND.EPS.GE.1D0) PRINT 7000,EPS IF(NP.EQ.-1.AND.EPS.LT.1D0) PRINT 7001,EPS IF(NP.GE.0) PRINT 7100,NP,EPS IF(NP.EQ.-2.AND.EPS.GE.1D0) PRINT 7150,EPS IF(NP.EQ.-2.AND.EPS.LT.1D0) PRINT 7151,EPS PRINT 7400, LAM,MRR,MRI write(6,7400) LAM,MRR,MRI PRINT 7200,DDELT write(6,7200) DDELT 7000 FORMAT('RANDOMLY ORIENTED OBLATE SPHEROIDS, A/B=',F11.7) 7001 FORMAT('RANDOMLY ORIENTED PROLATE SPHEROIDS, A/B=',F11.7) 7100 FORMAT('RANDOMLY ORIENTED CHEBYSHEV PARTICLES, T', & I1,'(',F5.2,')') 7150 FORMAT('RANDOMLY ORIENTED OBLATE CYLINDERS, D/L=',F11.7) 7151 FORMAT('RANDOMLY ORIENTED PROLATE CYLINDERS, D/L=',F11.7) 7200 FORMAT ('ACCURACY OF COMPUTATIONS DDELT = ',d8.2) 7400 FORMAT('LAM=',F10.6,3X,'MRR=',D10.4,3X,'MRI=',D10.4) DDELT=0.1D0*DDELT DO 600 IAX=1,NPNAX AXI=AXMAX-DAX*DFLOAT(IAX-1) R1=0.9999999d0*AXI R2=1.0000001D0*AXI NK=IDINT(AXI*NKMAX/AXMAX+2) IF (NK.GT.1000) PRINT 8001,NK IF (NK.GT.1000) STOP IF (NDISTR.EQ.3) CALL POWER (AXI,B,R1,R2) 8001 FORMAT ('NK=',I4,' I.E., IS GREATER THAN 1000. ', & 'EXECUTION TERMINATED.') CALL GAUSS (NK,0,0,XG,WG) Z1=(R2-R1)*0.5D0 Z2=(R1+R2)*0.5D0 Z3=R1*0.5D0 IF (NDISTR.EQ.5) GO TO 3 DO I=1,NK XG1(I)=Z1*XG(I)+Z2 WG1(I)=WG(I)*Z1 ENDDO GO TO 4 3 DO I=1,NK XG1(I)=Z3*XG(I)+Z3 WG1(I)=WG(I)*Z3 ENDDO DO I=NK+1,2*NK II=I-NK XG1(I)=Z1*XG(II)+Z2 WG1(I)=WG(II)*Z1 ENDDO NK=NK*2 4 CALL DISTRB (NK,XG1,WG1,NDISTR,AXI,B,GAM,R1,R2, & REFF,VEFF,PIN) PRINT 8002,R1,R2 8002 FORMAT('R1=',F10.6,' R2=',F10.6) IF (DABS(RAT-1D0).LE.1D-6) PRINT 8003, REFF,VEFF IF (DABS(RAT-1D0).GT.1D-6) PRINT 8004, REFF,VEFF 8003 FORMAT('EQUAL-VOLUME-SPHERE REFF=',F8.4,' VEFF=',F7.4) 8004 FORMAT('EQUAL-SURFACE-AREA-SPHERE REFF=',F8.4, & ' VEFF=',F7.4) PRINT 7250,NK write(6,7250) NK 7250 FORMAT('NUMBER OF GAUSSIAN QUADRATURE POINTS ', & 'IN SIZE AVERAGING =',I4) DO I=1,NPL ALPH1(I)=0D0 ALPH2(I)=0D0 ALPH3(I)=0D0 ALPH4(I)=0D0 BET1(I)=0D0 BET2(I)=0D0 ENDDO CSCAT=0D0 CEXTIN=0D0 L1MAX=0 DO 500 INK=1,NK I=NK-INK+1 A=RAT*XG1(I) XEV=2D0*PIN*A/DLAM IXXX=XEV+4.05D0*XEV**0.333333D0 INM1=MAX0(4,IXXX) IF (INM1.GE.NPN1) PRINT 7333, NPN1 IF (INM1.GE.NPN1) STOP 7333 FORMAT('CONVERGENCE IS NOT OBTAINED FOR NPN1=',I3, & '. EXECUTION TERMINATED') QEXT1=0D0 QSCA1=0D0 DO 50 NMA=INM1,NPN1 NMAX=NMA MMAX=1 NGAUSS=NMAX*NDGS IF (NGAUSS.GT.NPNG1) PRINT 7340, NGAUSS IF (NGAUSS.GT.NPNG1) STOP 7340 FORMAT('NGAUSS =',I3,' I.E. IS GREATER THAN NPNG1.', & ' EXECUTION TERMINATED') 7334 FORMAT(' NMAX =', I3,' DC2=',D8.2,' DC1=',D8.2) 7335 FORMAT(' NMAX1 =', I3,' DC2=',D8.2, & ' DC1=',D8.2) CALL CONST(NGAUSS,NMAX,MMAX,P,X,W,AN,ANN,S,SS,NP,EPS) CALL VARY(LAM,MRR,MRI,A,EPS,NP,NGAUSS,X,P,PPI,PIR,PII,R, & DR,DDR,DRR,DRI,NMAX) CALL TMATR0 (NGAUSS,X,W,AN,ANN,S,SS,PPI,PIR,PII,R,DR, & DDR,DRR,DRI,NMAX,NCHECK) QEXT=0D0 QSCA=0D0 DO N=1,NMAX N1=N+NMAX TR1NN=TR1(N,N) TI1NN=TI1(N,N) TR1NN1=TR1(N1,N1) TI1NN1=TI1(N1,N1) DN1=DFLOAT(2*N+1) QSCA=QSCA+DN1*(TR1NN*TR1NN+TI1NN*TI1NN & +TR1NN1*TR1NN1+TI1NN1*TI1NN1) QEXT=QEXT+(TR1NN+TR1NN1)*DN1 ENDDO DSCA=DABS((QSCA1-QSCA)/QSCA) DEXT=DABS((QEXT1-QEXT)/QEXT) C PRINT 7334, NMAX,DSCA,DEXT QEXT1=QEXT QSCA1=QSCA NMIN=DFLOAT(NMAX)/2D0+1D0 DO 10 N=NMIN,NMAX N1=N+NMAX TR1NN=TR1(N,N) TI1NN=TI1(N,N) TR1NN1=TR1(N1,N1) TI1NN1=TI1(N1,N1) DN1=DFLOAT(2*N+1) DQSCA=DN1*(TR1NN*TR1NN+TI1NN*TI1NN & +TR1NN1*TR1NN1+TI1NN1*TI1NN1) DQEXT=(TR1NN+TR1NN1)*DN1 DQSCA=DABS(DQSCA/QSCA) DQEXT=DABS(DQEXT/QEXT) NMAX1=N IF (DQSCA.LE.DDELT.AND.DQEXT.LE.DDELT) GO TO 12 10 CONTINUE 12 CONTINUE c PRINT 7335, NMAX1,DQSCA,DQEXT IF(DSCA.LE.DDELT.AND.DEXT.LE.DDELT) GO TO 55 IF (NMA.EQ.NPN1) PRINT 7333, NPN1 IF (NMA.EQ.NPN1) STOP 50 CONTINUE 55 NNNGGG=NGAUSS+1 IF (NGAUSS.EQ.NPNG1) PRINT 7336 MMAX=NMAX1 DO 150 NGAUS=NNNGGG,NPNG1 NGAUSS=NGAUS NGGG=2*NGAUSS 7336 FORMAT('WARNING: NGAUSS=NPNG1') 7337 FORMAT(' NG=',I3,' DC2=',D8.2,' DC1=',D8.2) CALL CONST(NGAUSS,NMAX,MMAX,P,X,W,AN,ANN,S,SS,NP,EPS) CALL VARY(LAM,MRR,MRI,A,EPS,NP,NGAUSS,X,P,PPI,PIR,PII,R, & DR,DDR,DRR,DRI,NMAX) CALL TMATR0 (NGAUSS,X,W,AN,ANN,S,SS,PPI,PIR,PII,R,DR, & DDR,DRR,DRI,NMAX,NCHECK) QEXT=0D0 QSCA=0D0 DO 104 N=1,NMAX N1=N+NMAX TR1NN=TR1(N,N) TI1NN=TI1(N,N) TR1NN1=TR1(N1,N1) TI1NN1=TI1(N1,N1) DN1=DFLOAT(2*N+1) QSCA=QSCA+DN1*(TR1NN*TR1NN+TI1NN*TI1NN & +TR1NN1*TR1NN1+TI1NN1*TI1NN1) QEXT=QEXT+(TR1NN+TR1NN1)*DN1 104 CONTINUE DSCA=DABS((QSCA1-QSCA)/QSCA) DEXT=DABS((QEXT1-QEXT)/QEXT) c PRINT 7337, NGGG,DSCA,DEXT QEXT1=QEXT QSCA1=QSCA IF(DSCA.LE.DDELT.AND.DEXT.LE.DDELT) GO TO 155 IF (NGAUS.EQ.NPNG1) PRINT 7336 150 CONTINUE 155 CONTINUE QSCA=0D0 QEXT=0D0 NNM=NMAX*2 DO 204 N=1,NNM QEXT=QEXT+TR1(N,N) 204 CONTINUE IF (NMAX1.GT.NPN4) PRINT 7550, NMAX1 7550 FORMAT ('nmax1 = ',I3, ', i.e. greater than NPN4.', & ' Execution terminated') IF (NMAX1.GT.NPN4) STOP DO 213 N2=1,NMAX1 NN2=N2+NMAX DO 213 N1=1,NMAX1 NN1=N1+NMAX ZZ1=TR1(N1,N2) RT11(1,N1,N2)=ZZ1 ZZ2=TI1(N1,N2) IT11(1,N1,N2)=ZZ2 ZZ3=TR1(N1,NN2) RT12(1,N1,N2)=ZZ3 ZZ4=TI1(N1,NN2) IT12(1,N1,N2)=ZZ4 ZZ5=TR1(NN1,N2) RT21(1,N1,N2)=ZZ5 ZZ6=TI1(NN1,N2) IT21(1,N1,N2)=ZZ6 ZZ7=TR1(NN1,NN2) RT22(1,N1,N2)=ZZ7 ZZ8=TI1(NN1,NN2) IT22(1,N1,N2)=ZZ8 QSCA=QSCA+ZZ1*ZZ1+ZZ2*ZZ2+ZZ3*ZZ3+ZZ4*ZZ4 & +ZZ5*ZZ5+ZZ6*ZZ6+ZZ7*ZZ7+ZZ8*ZZ8 213 CONTINUE C PRINT 7800,0,DABS(QEXT),QSCA,NMAX DO 220 M=1,NMAX1 CALL TMATR(M,NGAUSS,X,W,AN,ANN,S,SS,PPI,PIR,PII,R,DR, & DDR,DRR,DRI,NMAX,NCHECK) NM=NMAX-M+1 NM1=NMAX1-M+1 M1=M+1 QSC=0D0 DO 214 N2=1,NM1 NN2=N2+M-1 N22=N2+NM DO 214 N1=1,NM1 NN1=N1+M-1 N11=N1+NM ZZ1=TR1(N1,N2) RT11(M1,NN1,NN2)=ZZ1 ZZ2=TI1(N1,N2) IT11(M1,NN1,NN2)=ZZ2 ZZ3=TR1(N1,N22) RT12(M1,NN1,NN2)=ZZ3 ZZ4=TI1(N1,N22) IT12(M1,NN1,NN2)=ZZ4 ZZ5=TR1(N11,N2) RT21(M1,NN1,NN2)=ZZ5 ZZ6=TI1(N11,N2) IT21(M1,NN1,NN2)=ZZ6 ZZ7=TR1(N11,N22) RT22(M1,NN1,NN2)=ZZ7 ZZ8=TI1(N11,N22) IT22(M1,NN1,NN2)=ZZ8 QSC=QSC+(ZZ1*ZZ1+ZZ2*ZZ2+ZZ3*ZZ3+ZZ4*ZZ4 & +ZZ5*ZZ5+ZZ6*ZZ6+ZZ7*ZZ7+ZZ8*ZZ8)*2D0 214 CONTINUE NNM=2*NM QXT=0D0 c DO 215 N=1,NNM c QXT=QXT+TR1(N,N)*2D0 c 215 CONTINUE c QSCA=QSCA+QSC c QEXT=QEXT+QXT C PRINT 7800,M,DABS(QXT),QSC,NMAX 7800 FORMAT(' m=',I3,' qxt=',D12.6,' qsc=',D12.6, & ' nmax=',I3) 220 CONTINUE c COEFF1=LAM*LAM*0.5D0/P c CSCA=QSCA*COEFF1 c CEXT=-QEXT*COEFF1 c PRINT 7880, NMAX,NMAX1 7880 FORMAT ('nmax=',I3,' nmax1=',I3) c CALL GSP (NMAX1,CSCA,LAM,AL1,AL2,AL3,AL4,BE1,BE2,LMAX) c L1M=LMAX+1 c L1MAX=MAX(L1MAX,L1M) c WGII=WG1(I) c WGI=WGII*CSCA c DO 250 L1=1,L1M c ALPH1(L1)=ALPH1(L1)+AL1(L1)*WGI c ALPH2(L1)=ALPH2(L1)+AL2(L1)*WGI c ALPH3(L1)=ALPH3(L1)+AL3(L1)*WGI c ALPH4(L1)=ALPH4(L1)+AL4(L1)*WGI c BET1(L1)=BET1(L1)+BE1(L1)*WGI c BET2(L1)=BET2(L1)+BE2(L1)*WGI c 250 CONTINUE c CSCAT=CSCAT+WGI c CEXTIN=CEXTIN+CEXT*WGII C PRINT 6070, I,NMAX,NMAX1,NGAUSS 6070 FORMAT(4I6) 500 CONTINUE c DO 510 L1=1,L1MAX c ALPH1(L1)=ALPH1(L1)/CSCAT c ALPH2(L1)=ALPH2(L1)/CSCAT c ALPH3(L1)=ALPH3(L1)/CSCAT c ALPH4(L1)=ALPH4(L1)/CSCAT c BET1(L1)=BET1(L1)/CSCAT c BET2(L1)=BET2(L1)/CSCAT c 510 CONTINUE c WALB=CSCAT/CEXTIN c CALL HOVENR(L1MAX,ALPH1,ALPH2,ALPH3,ALPH4,BET1,BET2) c ASYMM=ALPH1(2)/3D0 c PRINT 9100,CEXTIN,CSCAT,WALB,ASYMM 9100 FORMAT('CEXT=',D12.6,2X,'CSCA=',D12.6,2X, & 2X,'W=',D12.6,2X,'=',D12.6) c write(6,9100) CEXTIN,CSCAT,WALB,ASYMM c temp=p*R1*R1 9101 FORMAT('QEXT=',D12.6,2X,'QSCA=',D12.6) c PRINT 9101,CEXTIN/temp,CSCAT/temp c write(6,9101) CEXTIN/temp,CSCAT/temp c IF (WALB.GT.1D0) PRINT 9111 9111 FORMAT ('WARNING: W IS GREATER THAN 1') c WRITE (10,580) WALB,L1MAX c DO L=1,L1MAX c WRITE (10,575) ALPH1(L),ALPH2(L),ALPH3(L),ALPH4(L), c & BET1(L),BET2(L) c ENDDO 575 FORMAT(6D14.7) 580 FORMAT(D14.8,I8) LMAX=L1MAX-1 c CALL MATR (ALPH1,ALPH2,ALPH3,ALPH4,BET1,BET2,LMAX,NPNA) 600 CONTINUE return END C********************************************************************** C * C INPUT PARAMETERS: * C * C NG = 2*NGAUSS - number of quadrature points on the * C interval (-1,1). NGAUSS.LE.NPNG1 * C NMAX,MMAX - maximum dimensions of the arrays. NMAX.LE.NPN1 * C MMAX.LE.NPN1 * C P - pi * C * C OUTPUT PARAMETERS: * C * C X,W - points and weights of the quadrature formula * C AN(N) = n*(n+1) * C ANN(N1,N2) = (1/2)*sqrt((2*n1+1)*(2*n2+1)/(n1*(n1+1)*n2*(n2+1))) * C S(I)=1/sin(arccos(x(i))) * C SS(I)=S(I)**2 * C * C********************************************************************** SUBROUTINE CONST (NGAUSS,NMAX,MMAX,P,X,W,AN,ANN,S,SS,NP,EPS) IMPLICIT REAL*16 (A-H,O-Z) INCLUDE 'tmq.par.f' REAL*16 X(NPNG2),W(NPNG2),X1(NPNG1),W1(NPNG1), * X2(NPNG1),W2(NPNG1), * S(NPNG2),SS(NPNG2), * AN(NPN1),ANN(NPN1,NPN1),DD(NPN1) DO 10 N=1,NMAX NN=N*(N+1) AN(N)=QFLOAT(NN) D=QSQRT(QFLOAT(2*N+1)/QFLOAT(NN)) DD(N)=D DO 10 N1=1,N DDD=D*DD(N1)*0.5Q0 ANN(N,N1)=DDD ANN(N1,N)=DDD 10 CONTINUE NG=2*NGAUSS IF (NP.EQ.-2) GO TO 11 CALL QGAUSS(NG,0,0,X,W) GO TO 19 11 NG1=DFLOAT(NGAUSS)/2D0 NG2=NGAUSS-NG1 XX=-QCOS(QATAN(EPS)) CALL QGAUSS(NG1,0,0,X1,W1) CALL QGAUSS(NG2,0,0,X2,W2) DO 12 I=1,NG1 W(I)=0.5Q0*(XX+1Q0)*W1(I) X(I)=0.5Q0*(XX+1Q0)*X1(I)+0.5Q0*(XX-1Q0) 12 CONTINUE DO 14 I=1,NG2 W(I+NG1)=-0.5Q0*XX*W2(I) X(I+NG1)=-0.5Q0*XX*X2(I)+0.5Q0*XX 14 CONTINUE DO 16 I=1,NGAUSS W(NG-I+1)=W(I) X(NG-I+1)=-X(I) 16 CONTINUE 19 DO 20 I=1,NGAUSS Y=X(I) Y=1Q0/(1Q0-Y*Y) SS(I)=Y SS(NG-I+1)=Y Y=QSQRT(Y) S(I)=Y S(NG-I+1)=Y 20 CONTINUE RETURN END C*************************************************************** SUBROUTINE QGAUSS ( N,IND1,IND2,Z,W ) IMPLICIT REAL*16 (A-H,P-Z) REAL*16 Z(N),W(N) DATA A,B,C /1Q0,2Q0,3Q0/ IND=MOD(N,2) K=N/2+IND F=QFLOAT(N) DO 100 I=1,K M=N+1-I IF(I.EQ.1) X=A-B/((F+A)*F) IF(I.EQ.2) X=(Z(N)-A)*4Q0+Z(N) IF(I.EQ.3) X=(Z(N-1)-Z(N))*1.6Q0+Z(N-1) IF(I.GT.3) X=(Z(M+1)-Z(M+2))*C+Z(M+3) IF(I.EQ.K.AND.IND.EQ.1) X=0Q0 NITER=0 CHECK=1Q-32 10 PB=1Q0 NITER=NITER+1 IF (NITER.LE.100) GO TO 15 c PRINT 5000, CHECK CHECK=CHECK*10Q0 15 PC=X DJ=A DO 20 J=2,N DJ=DJ+A PA=PB PB=PC 20 PC=X*PB+(X*PB-PA)*(DJ-A)/DJ PA=A/((PB-X*PC)*F) PB=PA*PC*(A-X*X) X=X-PB IF(QABS(PB).GT.CHECK*QABS(X)) GO TO 10 Z(M)=X W(M)=PA*PA*(A-X*X) IF(IND1.EQ.0) W(M)=B*W(M) IF(I.EQ.K.AND.IND.EQ.1) GO TO 100 Z(I)=-Z(M) W(I)=W(M) 100 CONTINUE 5000 format ('QGAUSS DOES NOT CONVERGE, CHECK=',f10.3) IF(IND2.NE.1) GO TO 110 PRINT 1100,N 1100 FORMAT(' *** POINTS AND WEIGHTS OF GAUSSIAN QUADRATURE FORMULA', * ' OF ',I4,'-TH ORDER') DO 105 I=1,K ZZ=-Z(I) 105 PRINT 1200,I,ZZ,I,W(I) 1200 FORMAT(' ',4X,'X(',I4,') = ',F17.14,5X,'W(',I4,') = ',F17.14) GO TO 115 110 CONTINUE C PRINT 1300,N 1300 FORMAT(' GAUSSIAN QUADRATURE FORMULA OF ',I4,'-TH ORDER IS USED') 115 CONTINUE IF(IND1.EQ.0) GO TO 140 DO 120 I=1,N 120 Z(I)=(A+Z(I))/B 140 CONTINUE RETURN END C********************************************************************** C * C INPUT PARAMETERS: * C * C LAM - wavelength of light * C MRR and MRI - real and imaginary parts of the refractive index * C A,EPS,NP - specify shape of the particle * C (see subroutines RSP1, RSP2, and RSP3) * C NG = NGAUSS*2 - number of gaussian quadrature points on the * C interval (-1,1) * C X - gaussian division points * C P - pi * C * C OUTPUT INFORMATION: * C * C PPI = PI**2 , where PI = (2*P)/LAM (wavenumber) * C PIR = PPI*MRR * C PII = PPI*MRI * C R and DR - see subroutines RSP1, RSP2, and RSP3 * C DDR=1/(PI*SQRT(R)) * C DRR+I*DRI=DDR/(MRR+I*MRI) * C NMAX - dimension of T(m)-matrices * C arrays J,Y,JR,JI,DJ,DY,DJR,DJI are transferred through * C COMMON /CBESS/ - see subroutine BESS * C * C********************************************************************** SUBROUTINE VARY (LAM,MRR,MRI,A,EPS,NP,NGAUSS,X,P,PPI,PIR,PII, * R,DR,DDR,DRR,DRI,NMAX) INCLUDE 'tmq.par.f' IMPLICIT REAL*16 (A-H,O-Z) REAL*16 X(NPNG2),R(NPNG2),DR(NPNG2),MRR,MRI,LAM, * Z(NPNG2),ZR(NPNG2),ZI(NPNG2), * J(NPNG2,NPN1),Y(NPNG2,NPN1),JR(NPNG2,NPN1), * JI(NPNG2,NPN1),DJ(NPNG2,NPN1), * DJR(NPNG2,NPN1),DJI(NPNG2,NPN1),DDR(NPNG2), * DRR(NPNG2),DRI(NPNG2), * DY(NPNG2,NPN1) COMMON /CBESS/ J,Y,JR,JI,DJ,DY,DJR,DJI NG=NGAUSS*2 IF (NP.EQ.-1) CALL RSP1(X,NG,NGAUSS,A,EPS,NP,R,DR) IF (NP.GE.0) CALL RSP2(X,NG,A,EPS,NP,R,DR) IF (NP.EQ.-2) CALL RSP3(X,NG,NGAUSS,A,EPS,R,DR) PI=P*2Q0/LAM PPI=PI*PI PIR=PPI*MRR PII=PPI*MRI V=1Q0/(MRR*MRR+MRI*MRI) PRR=MRR*V PRI=-MRI*V TA=0Q0 DO 10 I=1,NG VV=QSQRT(R(I)) V=VV*PI TA=MAX(TA,V) VV=1Q0/V DDR(I)=VV DRR(I)=PRR*VV DRI(I)=PRI*VV V1=V*MRR V2=V*MRI Z(I)=V ZR(I)=V1 ZI(I)=V2 10 CONTINUE IF (NMAX.GT.NPN1) PRINT 9000,NMAX,NPN1 IF (NMAX.GT.NPN1) STOP 9000 FORMAT(' NMAX = ',I2,', i.e., greater than ',I3) TB=TA*QSQRT(MRR*MRR+MRI*MRI) TB=QMAX1(TB,QFLOAT(NMAX)) NNMAX1=8.0Q0*QSQRT(QMAX1(TA,QFLOAT(NMAX)))+3Q0 NNMAX2=(TB+4Q0*(TB**0.33333Q0)+8.0Q0*QSQRT(TB)) NNMAX2=NNMAX2-NMAX+5 CALL BESS(Z,ZR,ZI,NG,NMAX,NNMAX1,NNMAX2) RETURN END C********************************************************************** C * C Calculation of the functions R(I)=r(y(I))**2 and * C DR(I)=((d/dy)r(y))/r(y) and horizontal semi-axis A * C for a spheroid specified by the parameters REV (equal-volume- * C sphere radius) and EPS=A/B (ratio of the semi-axes). * C Y(I)=arccos(X(I)) * C 1.LE.I.LE.NG * C X - arguments * C * C********************************************************************** SUBROUTINE RSP1 (X,NG,NGAUSS,REV,EPS,NP,R,DR) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 X(NG),R(NG),DR(NG) A=REV*EPS**(1Q0/3Q0) AA=A*A EE=EPS*EPS EE1=EE-1Q0 DO 50 I=1,NGAUSS C=X(I) CC=C*C SS=1Q0-CC S=QSQRT(SS) RR=1Q0/(SS+EE*CC) R(I)=AA*RR R(NG-I+1)=R(I) DR(I)=RR*C*S*EE1 DR(NG-I+1)=-DR(I) 50 CONTINUE RETURN END C********************************************************************** C * C Calculation of the functions R(I)=r(y(i))**2 and * C DR(I)=((d/dy)r(y))/r(y) and parameter R0 for a Chebyshev * C particle specified by the parameters REV (equal-volume-sphere * C radius), EPS, and N. * C Y(I)=arccos(X(I)) * C 1.LE.I.LE.NG * C X - arguments * C * C********************************************************************** SUBROUTINE RSP2 (X,NG,REV,EPS,N,R,DR) IMPLICIT REAL*16 (A-H,O-Z) double precision temp0 REAL*16 X(NG),R(NG),DR(NG) DNP=QFLOAT(N) DN=DNP*DNP DN4=DN*4Q0 EP=EPS*EPS A=1Q0+1.5Q0*EP*(DN4-2Q0)/(DN4-1Q0) I=(DNP+0.1Q0)*0.5Q0 I=2*I IF (I.EQ.N) A=A-3Q0*EPS*(1Q0+0.25Q0*EP)/ * (DN-1Q0)-0.25Q0*EP*EPS/(9Q0*DN-1Q0) R0=REV*A**(-1Q0/3Q0) DO 50 I=1,NG temp0=X(I) temp=dacos(temp0) XI=temp*DNP RI=R0*(1Q0+EPS*QCOS(XI)) R(I)=RI*RI DR(I)=-R0*EPS*DNP*QSIN(XI)/RI 50 CONTINUE RETURN END C********************************************************************** C * C Calculation of the functions R(I)=R(Y(I))**2 and * C DR(I)=((d/dy)r(y))/r(y) * C for a cylinder specified by the parameters REV (equal-volume- * C sphere radius) and EPS=A/H (ratio of radius to semi-length) * C Y(I)=arccos(X(I)) * C 1.LE.I.LE.NG * C X - arguments * C * C********************************************************************** SUBROUTINE RSP3 (X,NG,NGAUSS,REV,EPS,R,DR) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 X(NG),R(NG),DR(NG) H=REV*( (2Q0/(3Q0*EPS*EPS))**(1Q0/3Q0) ) A=H*EPS DO 50 I=1,NGAUSS CO=-X(I) SI=QSQRT(1Q0-CO*CO) IF (SI/CO.GT.A/H) GO TO 20 RAD=H/CO RTHET=H*SI/(CO*CO) GO TO 30 20 RAD=A/SI RTHET=-A*CO/(SI*SI) 30 R(I)=RAD*RAD R(NG-I+1)=R(I) DR(I)=-RTHET/RAD DR(NG-I+1)=-DR(I) 50 CONTINUE RETURN END C********************************************************************** C * C Calculation of spherical Bessel functions of the first kind * C J(I,N) = j_n(x) and second kind Y(I,N) = y_n(x) * C of real-valued argument X(I) and first kind JR(I,N)+I*JI(I,N) = * C = j_n(z) of complex argument Z(I)=XR(I)+I*XI(I), as well as * C the functions * C * C DJ(I,N) = (1/x)(d/dx)(x*j_n(x)) , * C DY(I,N) = (1/x)(d/dx)(x*y_n(x)) , * C DJR(I,N) = Re ((1/z)(d/dz)(z*j_n(x)) , * C DJI(I,N) = Im ((1/z)(d/dz)(z*j_n(x)) . * C * C 1.LE.N.LE.NMAX * C NMAX.LE.NPN1 * C X,XR,XI - arguments * C 1.LE.I.LE.NG * C Arrays J,Y,JR,JI,DJ,DY,DJR,DJI are in * C COMMON /CBESS/ * C Parameters NNMAX1 and NMAX2 determine computational accuracy * C * C********************************************************************** SUBROUTINE BESS (X,XR,XI,NG,NMAX,NNMAX1,NNMAX2) INCLUDE 'tmq.par.f' IMPLICIT REAL*16 (A-H,O-Z) REAL*16 X(NG),XR(NG),XI(NG), * J(NPNG2,NPN1),Y(NPNG2,NPN1),JR(NPNG2,NPN1), * JI(NPNG2,NPN1),DJ(NPNG2,NPN1),DY(NPNG2,NPN1), * DJR(NPNG2,NPN1),DJI(NPNG2,NPN1), * AJ(NPN1),AY(NPN1),AJR(NPN1),AJI(NPN1), * ADJ(NPN1),ADY(NPN1),ADJR(NPN1), * ADJI(NPN1) COMMON /CBESS/ J,Y,JR,JI,DJ,DY,DJR,DJI DO 10 I=1,NG XX=X(I) CALL RJB(XX,AJ,ADJ,NMAX,NNMAX1) CALL RYB(XX,AY,ADY,NMAX) YR=XR(I) YI=XI(I) CALL CJB(YR,YI,AJR,AJI,ADJR,ADJI,NMAX,NNMAX2) DO 10 N=1,NMAX J(I,N)=AJ(N) Y(I,N)=AY(N) JR(I,N)=AJR(N) JI(I,N)=AJI(N) DJ(I,N)=ADJ(N) DY(I,N)=ADY(N) DJR(I,N)=ADJR(N) DJI(I,N)=ADJI(N) 10 CONTINUE RETURN END C********************************************************************** C * C Calculation of spherical Bessel functions of the first kind j * C of real-valued argument x of orders from 1 to NMAX by using * C backward recursion. Parametr NNMAX determines numerical accuracy. * C U - function (1/x)(d/dx)(x*j(x)) * C * C********************************************************************** SUBROUTINE RJB(X,Y,U,NMAX,NNMAX) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 Y(NMAX),U(NMAX),Z(900) L=NMAX+NNMAX XX=1Q0/X Z(L)=1Q0/(QFLOAT(2*L+1)*XX) L1=L-1 DO 5 I=1,L1 I1=L-I Z(I1)=1Q0/(QFLOAT(2*I1+1)*XX-Z(I1+1)) 5 CONTINUE Z0=1Q0/(XX-Z(1)) Y0=Z0*QCOS(X)*XX Y1=Y0*Z(1) U(1)=Y0-Y1*XX Y(1)=Y1 DO 10 I=2,NMAX YI1=Y(I-1) YI=YI1*Z(I) U(I)=YI1-QFLOAT(I)*YI*XX Y(I)=YI 10 CONTINUE RETURN END C********************************************************************** C * C Calculation of spherical Bessel functions of the second kind y * C of real-valued argument x of orders from 1 to NMAX by using upward* C recursion. V - function (1/x)(d/dx)(x*y(x)) * C * C********************************************************************** SUBROUTINE RYB(X,Y,V,NMAX) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 Y(NMAX),V(NMAX) C=QCOS(X) S=QSIN(X) X1=1Q0/X X2=X1*X1 X3=X2*X1 Y1=-C*X2-S*X1 Y(1)=Y1 Y(2)=(-3Q0*X3+X1)*C-3Q0*X2*S NMAX1=NMAX-1 DO 5 I=2,NMAX1 5 Y(I+1)=QFLOAT(2*I+1)*X1*Y(I)-Y(I-1) V(1)=-X1*(C+Y1) DO 10 I=2,NMAX 10 V(I)=Y(I-1)-QFLOAT(I)*X1*Y(I) RETURN END C********************************************************************** C * C Calculation of spherical Bessel functions of the first kind * C j=JR+I*JI of complex argument x=XR+I*XI of orders from 1 to NMAX * C by using backward recursion. Parametr NNMAX determines numerical * C accuracy. U=UR+I*UI - function (1/x)(d/dx)(x*j(x)) * C * C********************************************************************** SUBROUTINE CJB (XR,XI,YR,YI,UR,UI,NMAX,NNMAX) INCLUDE 'tmq.par.f' IMPLICIT REAL*16 (A-H,O-Z) REAL*16 YR(NMAX),YI(NMAX),UR(NMAX),UI(NMAX) REAL*16 CYR(NPN1),CYI(NPN1),CZR(1200),CZI(1200), * CUR(NPN1),CUI(NPN1) L=NMAX+NNMAX XRXI=1Q0/(XR*XR+XI*XI) CXXR=XR*XRXI CXXI=-XI*XRXI QF=1Q0/QFLOAT(2*L+1) CZR(L)=XR*QF CZI(L)=XI*QF L1=L-1 DO 5 I=1,L1 I1=L-I QF=QFLOAT(2*I1+1) AR=QF*CXXR-CZR(I1+1) AI=QF*CXXI-CZI(I1+1) ARI=1Q0/(AR*AR+AI*AI) CZR(I1)=AR*ARI CZI(I1)=-AI*ARI 5 CONTINUE AR=CXXR-CZR(1) AI=CXXI-CZI(1) ARI=1Q0/(AR*AR+AI*AI) CZ0R=AR*ARI CZ0I=-AI*ARI CR=QCOS(XR)*QCOSH(XI) CI=-QSIN(XR)*QSINH(XI) AR=CZ0R*CR-CZ0I*CI AI=CZ0I*CR+CZ0R*CI CY0R=AR*CXXR-AI*CXXI CY0I=AI*CXXR+AR*CXXI CY1R=CY0R*CZR(1)-CY0I*CZI(1) CY1I=CY0I*CZR(1)+CY0R*CZI(1) AR=CY1R*CXXR-CY1I*CXXI AI=CY1I*CXXR+CY1R*CXXI CU1R=CY0R-AR CU1I=CY0I-AI CYR(1)=CY1R CYI(1)=CY1I CUR(1)=CU1R CUI(1)=CU1I YR(1)=CY1R YI(1)=CY1I UR(1)=CU1R UI(1)=CU1I DO 10 I=2,NMAX QI=QFLOAT(I) CYI1R=CYR(I-1) CYI1I=CYI(I-1) CYIR=CYI1R*CZR(I)-CYI1I*CZI(I) CYII=CYI1I*CZR(I)+CYI1R*CZI(I) AR=CYIR*CXXR-CYII*CXXI AI=CYII*CXXR+CYIR*CXXI CUIR=CYI1R-QI*AR CUII=CYI1I-QI*AI CYR(I)=CYIR CYI(I)=CYII CUR(I)=CUIR CUI(I)=CUII YR(I)=CYIR YI(I)=CYII UR(I)=CUIR UI(I)=CUII 10 CONTINUE RETURN END C********************************************************************** C * C calculation of the T(0) matrix for an axially symmetric particle * C * C Output information: * C * C Arrays TR1 and TI1 (real and imaginary parts of the * C T(0) matrix) are in COMMON /CT/ * C * C********************************************************************** SUBROUTINE TMATR0 (NGAUSS,X,W,AN,ANN,S,SS,PPI,PIR,PII,R,DR,DDR, * DRR,DRI,NMAX,NCHECK) INCLUDE 'tmq.par.f' IMPLICIT REAL*16 (A-H,O-Z) REAL*16 X(NPNG2),W(NPNG2),AN(NPN1),S(NPNG2),SS(NPNG2), * R(NPNG2),DR(NPNG2),SIG(NPN2), * J(NPNG2,NPN1),Y(NPNG2,NPN1), * JR(NPNG2,NPN1),JI(NPNG2,NPN1),DJ(NPNG2,NPN1), * DY(NPNG2,NPN1),DJR(NPNG2,NPN1), * DJI(NPNG2,NPN1),DDR(NPNG2),DRR(NPNG2), * D1(NPNG2,NPN1),D2(NPNG2,NPN1), * DRI(NPNG2),DS(NPNG2),DSS(NPNG2),RR(NPNG2), * DV1(NPN1),DV2(NPN1) REAL*16 R11(NPN1,NPN1),R12(NPN1,NPN1), * R21(NPN1,NPN1),R22(NPN1,NPN1), * I11(NPN1,NPN1),I12(NPN1,NPN1), * I21(NPN1,NPN1),I22(NPN1,NPN1), * RG11(NPN1,NPN1),RG12(NPN1,NPN1), * RG21(NPN1,NPN1),RG22(NPN1,NPN1), * IG11(NPN1,NPN1),IG12(NPN1,NPN1), * IG21(NPN1,NPN1),IG22(NPN1,NPN1), * ANN(NPN1,NPN1), * QR(NPN2,NPN2),QI(NPN2,NPN2), * RGQR(NPN2,NPN2),RGQI(NPN2,NPN2), * TQR(NPN2,NPN2),TQI(NPN2,NPN2), * TRGQR(NPN2,NPN2),TRGQI(NPN2,NPN2) double precision TR1(NPN2,NPN2),TI1(NPN2,NPN2) double precision PLUS(NPN6*NPN4*NPN4*8) COMMON /TMAT/ PLUS, & R11,R12,R21,R22,I11,I12,I21,I22,RG11,RG12,RG21,RG22, & IG11,IG12,IG21,IG22 COMMON /CBESS/ J,Y,JR,JI,DJ,DY,DJR,DJI COMMON /CT/ TR1,TI1 COMMON /CTT/ QR,QI,RGQR,RGQI MM1=1 NNMAX=NMAX+NMAX NG=2*NGAUSS NGSS=NG FACTOR=1Q0 IF (NCHECK.EQ.1) THEN NGSS=NGAUSS FACTOR=2Q0 ELSE CONTINUE ENDIF SI=1Q0 DO 5 N=1,NNMAX SI=-SI SIG(N)=SI 5 CONTINUE 20 DO 25 I=1,NGAUSS I1=NGAUSS+I I2=NGAUSS-I+1 CALL VIG (X(I1),NMAX,0,DV1,DV2) DO 25 N=1,NMAX SI=SIG(N) DD1=DV1(N) DD2=DV2(N) D1(I1,N)=DD1 D2(I1,N)=DD2 D1(I2,N)=DD1*SI D2(I2,N)=-DD2*SI 25 CONTINUE 30 DO 40 I=1,NGSS RR(I)=W(I)*R(I) 40 CONTINUE DO 300 N1=MM1,NMAX AN1=AN(N1) DO 300 N2=MM1,NMAX AN2=AN(N2) AR12=0Q0 AR21=0Q0 AI12=0Q0 AI21=0Q0 GR12=0Q0 GR21=0Q0 GI12=0Q0 GI21=0Q0 IF (NCHECK.EQ.1.AND.SIG(N1+N2).LT.0Q0) GO TO 205 DO 200 I=1,NGSS D1N1=D1(I,N1) D2N1=D2(I,N1) D1N2=D1(I,N2) D2N2=D2(I,N2) A12=D1N1*D2N2 A21=D2N1*D1N2 A22=D2N1*D2N2 AA1=A12+A21 QJ1=J(I,N1) QY1=Y(I,N1) QJR2=JR(I,N2) QJI2=JI(I,N2) QDJR2=DJR(I,N2) QDJI2=DJI(I,N2) QDJ1=DJ(I,N1) QDY1=DY(I,N1) C1R=QJR2*QJ1 C1I=QJI2*QJ1 B1R=C1R-QJI2*QY1 B1I=C1I+QJR2*QY1 C2R=QJR2*QDJ1 C2I=QJI2*QDJ1 B2R=C2R-QJI2*QDY1 B2I=C2I+QJR2*QDY1 DDRI=DDR(I) C3R=DDRI*C1R C3I=DDRI*C1I B3R=DDRI*B1R B3I=DDRI*B1I C4R=QDJR2*QJ1 C4I=QDJI2*QJ1 B4R=C4R-QDJI2*QY1 B4I=C4I+QDJR2*QY1 DRRI=DRR(I) DRII=DRI(I) C5R=C1R*DRRI-C1I*DRII C5I=C1I*DRRI+C1R*DRII B5R=B1R*DRRI-B1I*DRII B5I=B1I*DRRI+B1R*DRII URI=DR(I) RRI=RR(I) F1=RRI*A22 F2=RRI*URI*AN1*A12 AR12=AR12+F1*B2R+F2*B3R AI12=AI12+F1*B2I+F2*B3I GR12=GR12+F1*C2R+F2*C3R GI12=GI12+F1*C2I+F2*C3I F2=RRI*URI*AN2*A21 AR21=AR21+F1*B4R+F2*B5R AI21=AI21+F1*B4I+F2*B5I GR21=GR21+F1*C4R+F2*C5R GI21=GI21+F1*C4I+F2*C5I 200 CONTINUE 205 AN12=ANN(N1,N2)*FACTOR R12(N1,N2)=AR12*AN12 R21(N1,N2)=AR21*AN12 I12(N1,N2)=AI12*AN12 I21(N1,N2)=AI21*AN12 RG12(N1,N2)=GR12*AN12 RG21(N1,N2)=GR21*AN12 IG12(N1,N2)=GI12*AN12 IG21(N1,N2)=GI21*AN12 300 CONTINUE TPIR=PIR TPII=PII TPPI=PPI NM=NMAX DO 310 N1=MM1,NMAX K1=N1-MM1+1 KK1=K1+NM DO 310 N2=MM1,NMAX K2=N2-MM1+1 KK2=K2+NM TAR12= I12(N1,N2) TAI12=-R12(N1,N2) TGR12= IG12(N1,N2) TGI12=-RG12(N1,N2) TAR21=-I21(N1,N2) TAI21= R21(N1,N2) TGR21=-IG21(N1,N2) TGI21= RG21(N1,N2) TQR(K1,K2)=TPIR*TAR21-TPII*TAI21+TPPI*TAR12 TQI(K1,K2)=TPIR*TAI21+TPII*TAR21+TPPI*TAI12 TRGQR(K1,K2)=TPIR*TGR21-TPII*TGI21+TPPI*TGR12 TRGQI(K1,K2)=TPIR*TGI21+TPII*TGR21+TPPI*TGI12 TQR(K1,KK2)=0Q0 TQI(K1,KK2)=0Q0 TRGQR(K1,KK2)=0Q0 TRGQI(K1,KK2)=0Q0 TQR(KK1,K2)=0Q0 TQI(KK1,K2)=0Q0 TRGQR(KK1,K2)=0Q0 TRGQI(KK1,K2)=0Q0 TQR(KK1,KK2)=TPIR*TAR12-TPII*TAI12+TPPI*TAR21 TQI(KK1,KK2)=TPIR*TAI12+TPII*TAR12+TPPI*TAI21 TRGQR(KK1,KK2)=TPIR*TGR12-TPII*TGI12+TPPI*TGR21 TRGQI(KK1,KK2)=TPIR*TGI12+TPII*TGR12+TPPI*TGI21 310 CONTINUE NNMAX=2*NM DO 320 N1=1,NNMAX DO 320 N2=1,NNMAX QR(N1,N2)=TQR(N1,N2) QI(N1,N2)=TQI(N1,N2) RGQR(N1,N2)=TRGQR(N1,N2) RGQI(N1,N2)=TRGQI(N1,N2) 320 CONTINUE CALL TT(NMAX,NCHECK) RETURN END C********************************************************************** C * C Calculation of the T(M) matrix, M.GE.1, for an axially symmetric * C particle * C * C Input parameters: * C * C M.GE.1 * C NG = NGAUSS*2 - number of gaussian division points on the * C interval (-1,1) * C W - quadrature weights * C AN,ANN - see subroutine CONST * C S,SS - see subroutine CONST * C ARRAYS DV1,DV2,DV3,DV4 are in COMMON /DV/ - * C see subroutine DVIG * C PPI,PIR,PII - see subroutine VARY * C R J DR - see subroutines RSP1 and RSP2 * C DDR,DRR,DRI - see subroutine VARY * C NMAX - dimension of the T(M) matrix * C Arrays J,Y,JR,JI,DJ,DY,DJR,DJI are in * C COMMON /CBESS/ - see subroutine BESS * C * C Output parameters: * C * C Arrays TR1,TI1 (real and imaginary parts of the T(M) matrix) * C are in COMMON /CT/ * C * C********************************************************************** SUBROUTINE TMATR (M,NGAUSS,X,W,AN,ANN,S,SS,PPI,PIR,PII,R,DR,DDR, * DRR,DRI,NMAX,NCHECK) INCLUDE 'tmq.par.f' IMPLICIT REAL*16 (A-H,O-Z) REAL*16 X(NPNG2),W(NPNG2),AN(NPN1),S(NPNG2),SS(NPNG2), * R(NPNG2),DR(NPNG2),SIG(NPN2), * J(NPNG2,NPN1),Y(NPNG2,NPN1), * JR(NPNG2,NPN1),JI(NPNG2,NPN1),DJ(NPNG2,NPN1), * DY(NPNG2,NPN1),DJR(NPNG2,NPN1), * DJI(NPNG2,NPN1),DDR(NPNG2),DRR(NPNG2), * D1(NPNG2,NPN1),D2(NPNG2,NPN1), * DRI(NPNG2),DS(NPNG2),DSS(NPNG2),RR(NPNG2), * DV1(NPN1),DV2(NPN1) REAL*16 R11(NPN1,NPN1),R12(NPN1,NPN1), * R21(NPN1,NPN1),R22(NPN1,NPN1), * I11(NPN1,NPN1),I12(NPN1,NPN1), * I21(NPN1,NPN1),I22(NPN1,NPN1), * RG11(NPN1,NPN1),RG12(NPN1,NPN1), * RG21(NPN1,NPN1),RG22(NPN1,NPN1), * IG11(NPN1,NPN1),IG12(NPN1,NPN1), * IG21(NPN1,NPN1),IG22(NPN1,NPN1), * ANN(NPN1,NPN1), * QR(NPN2,NPN2),QI(NPN2,NPN2), * RGQR(NPN2,NPN2),RGQI(NPN2,NPN2), * TQR(NPN2,NPN2),TQI(NPN2,NPN2), * TRGQR(NPN2,NPN2),TRGQI(NPN2,NPN2) double precision TR1(NPN2,NPN2),TI1(NPN2,NPN2) double precision PLUS(NPN6*NPN4*NPN4*8) COMMON /TMAT/ PLUS, & R11,R12,R21,R22,I11,I12,I21,I22,RG11,RG12,RG21,RG22, & IG11,IG12,IG21,IG22 COMMON /CBESS/ J,Y,JR,JI,DJ,DY,DJR,DJI COMMON /CT/ TR1,TI1 COMMON /CTT/ QR,QI,RGQR,RGQI MM1=M QM=QFLOAT(M) QMM=QM*QM NG=2*NGAUSS NGSS=NG FACTOR=1Q0 IF (NCHECK.EQ.1) THEN NGSS=NGAUSS FACTOR=2Q0 ELSE CONTINUE ENDIF SI=1Q0 NM=NMAX+NMAX DO 5 N=1,NM SI=-SI SIG(N)=SI 5 CONTINUE 20 DO 25 I=1,NGAUSS I1=NGAUSS+I I2=NGAUSS-I+1 CALL VIG (X(I1),NMAX,M,DV1,DV2) DO 25 N=1,NMAX SI=SIG(N+M) DD1=DV1(N) DD2=DV2(N) D1(I1,N)=DD1 D2(I1,N)=DD2 D1(I2,N)=DD1*SI D2(I2,N)=-DD2*SI 25 CONTINUE 30 DO 40 I=1,NGSS WR=W(I)*R(I) DS(I)=S(I)*QM*WR DSS(I)=SS(I)*QMM RR(I)=WR 40 CONTINUE DO 300 N1=MM1,NMAX AN1=AN(N1) DO 300 N2=MM1,NMAX AN2=AN(N2) AR11=0Q0 AR12=0Q0 AR21=0Q0 AR22=0Q0 AI11=0Q0 AI12=0Q0 AI21=0Q0 AI22=0Q0 GR11=0Q0 GR12=0Q0 GR21=0Q0 GR22=0Q0 GI11=0Q0 GI12=0Q0 GI21=0Q0 GI22=0Q0 SI=SIG(N1+N2) DO 200 I=1,NGSS D1N1=D1(I,N1) D2N1=D2(I,N1) D1N2=D1(I,N2) D2N2=D2(I,N2) A11=D1N1*D1N2 A12=D1N1*D2N2 A21=D2N1*D1N2 A22=D2N1*D2N2 AA1=A12+A21 AA2=A11*DSS(I)+A22 QJ1=J(I,N1) QY1=Y(I,N1) QJR2=JR(I,N2) QJI2=JI(I,N2) QDJR2=DJR(I,N2) QDJI2=DJI(I,N2) QDJ1=DJ(I,N1) QDY1=DY(I,N1) C1R=QJR2*QJ1 C1I=QJI2*QJ1 B1R=C1R-QJI2*QY1 B1I=C1I+QJR2*QY1 C2R=QJR2*QDJ1 C2I=QJI2*QDJ1 B2R=C2R-QJI2*QDY1 B2I=C2I+QJR2*QDY1 DDRI=DDR(I) C3R=DDRI*C1R C3I=DDRI*C1I B3R=DDRI*B1R B3I=DDRI*B1I C4R=QDJR2*QJ1 C4I=QDJI2*QJ1 B4R=C4R-QDJI2*QY1 B4I=C4I+QDJR2*QY1 DRRI=DRR(I) DRII=DRI(I) C5R=C1R*DRRI-C1I*DRII C5I=C1I*DRRI+C1R*DRII B5R=B1R*DRRI-B1I*DRII B5I=B1I*DRRI+B1R*DRII C6R=QDJR2*QDJ1 C6I=QDJI2*QDJ1 B6R=C6R-QDJI2*QDY1 B6I=C6I+QDJR2*QDY1 C7R=C4R*DDRI C7I=C4I*DDRI B7R=B4R*DDRI B7I=B4I*DDRI C8R=C2R*DRRI-C2I*DRII C8I=C2I*DRRI+C2R*DRII B8R=B2R*DRRI-B2I*DRII B8I=B2I*DRRI+B2R*DRII URI=DR(I) DSI=DS(I) DSSI=DSS(I) RRI=RR(I) IF (NCHECK.EQ.1.AND.SI.GT.0Q0) GO TO 150 E1=DSI*AA1 AR11=AR11+E1*B1R AI11=AI11+E1*B1I GR11=GR11+E1*C1R GI11=GI11+E1*C1I IF (NCHECK.EQ.1) GO TO 160 150 F1=RRI*AA2 F2=RRI*URI*AN1*A12 AR12=AR12+F1*B2R+F2*B3R AI12=AI12+F1*B2I+F2*B3I GR12=GR12+F1*C2R+F2*C3R GI12=GI12+F1*C2I+F2*C3I F2=RRI*URI*AN2*A21 AR21=AR21+F1*B4R+F2*B5R AI21=AI21+F1*B4I+F2*B5I GR21=GR21+F1*C4R+F2*C5R GI21=GI21+F1*C4I+F2*C5I IF (NCHECK.EQ.1) GO TO 200 160 E2=DSI*URI*A11 E3=E2*AN2 E2=E2*AN1 AR22=AR22+E1*B6R+E2*B7R+E3*B8R AI22=AI22+E1*B6I+E2*B7I+E3*B8I GR22=GR22+E1*C6R+E2*C7R+E3*C8R GI22=GI22+E1*C6I+E2*C7I+E3*C8I 200 CONTINUE AN12=ANN(N1,N2)*FACTOR R11(N1,N2)=AR11*AN12 R12(N1,N2)=AR12*AN12 R21(N1,N2)=AR21*AN12 R22(N1,N2)=AR22*AN12 I11(N1,N2)=AI11*AN12 I12(N1,N2)=AI12*AN12 I21(N1,N2)=AI21*AN12 I22(N1,N2)=AI22*AN12 RG11(N1,N2)=GR11*AN12 RG12(N1,N2)=GR12*AN12 RG21(N1,N2)=GR21*AN12 RG22(N1,N2)=GR22*AN12 IG11(N1,N2)=GI11*AN12 IG12(N1,N2)=GI12*AN12 IG21(N1,N2)=GI21*AN12 IG22(N1,N2)=GI22*AN12 300 CONTINUE TPIR=PIR TPII=PII TPPI=PPI NM=NMAX-MM1+1 DO 310 N1=MM1,NMAX K1=N1-MM1+1 KK1=K1+NM DO 310 N2=MM1,NMAX K2=N2-MM1+1 KK2=K2+NM TAR11=-R11(N1,N2) TAI11=-I11(N1,N2) TGR11=-RG11(N1,N2) TGI11=-IG11(N1,N2) TAR12= I12(N1,N2) TAI12=-R12(N1,N2) TGR12= IG12(N1,N2) TGI12=-RG12(N1,N2) TAR21=-I21(N1,N2) TAI21= R21(N1,N2) TGR21=-IG21(N1,N2) TGI21= RG21(N1,N2) TAR22=-R22(N1,N2) TAI22=-I22(N1,N2) TGR22=-RG22(N1,N2) TGI22=-IG22(N1,N2) TQR(K1,K2)=TPIR*TAR21-TPII*TAI21+TPPI*TAR12 TQI(K1,K2)=TPIR*TAI21+TPII*TAR21+TPPI*TAI12 TRGQR(K1,K2)=TPIR*TGR21-TPII*TGI21+TPPI*TGR12 TRGQI(K1,K2)=TPIR*TGI21+TPII*TGR21+TPPI*TGI12 TQR(K1,KK2)=TPIR*TAR11-TPII*TAI11+TPPI*TAR22 TQI(K1,KK2)=TPIR*TAI11+TPII*TAR11+TPPI*TAI22 TRGQR(K1,KK2)=TPIR*TGR11-TPII*TGI11+TPPI*TGR22 TRGQI(K1,KK2)=TPIR*TGI11+TPII*TGR11+TPPI*TGI22 TQR(KK1,K2)=TPIR*TAR22-TPII*TAI22+TPPI*TAR11 TQI(KK1,K2)=TPIR*TAI22+TPII*TAR22+TPPI*TAI11 TRGQR(KK1,K2)=TPIR*TGR22-TPII*TGI22+TPPI*TGR11 TRGQI(KK1,K2)=TPIR*TGI22+TPII*TGR22+TPPI*TGI11 TQR(KK1,KK2)=TPIR*TAR12-TPII*TAI12+TPPI*TAR21 TQI(KK1,KK2)=TPIR*TAI12+TPII*TAR12+TPPI*TAI21 TRGQR(KK1,KK2)=TPIR*TGR12-TPII*TGI12+TPPI*TGR21 TRGQI(KK1,KK2)=TPIR*TGI12+TPII*TGR12+TPPI*TGI21 310 CONTINUE NNMAX=2*NM DO 320 N1=1,NNMAX DO 320 N2=1,NNMAX QR(N1,N2)=TQR(N1,N2) QI(N1,N2)=TQI(N1,N2) RGQR(N1,N2)=TRGQR(N1,N2) RGQI(N1,N2)=TRGQI(N1,N2) 320 CONTINUE CALL TT(NM,NCHECK) RETURN END C***************************************************************** c c Calculation of the functiONS c DV1(n)=dvig(0,m,n,arccos x) c and c DV2(n)=[d/d(arccos x)] dvig(0,m,n,arccos x) c 1.LE.N.LE.NMAX c 0.LE.x.LE.1 SUBROUTINE VIG (X,NMAX,M,DV1,DV2) INCLUDE 'tmq.par.f' IMPLICIT REAL*16 (A-H,O-Z) REAL*16 DV1(NPN1), DV2(NPN1) A=1Q0 QS=QSQRT(1Q0-X*X) QS1=1Q0/QS DO 1 N=1,NMAX DV1(N)=0Q0 DV2(N)=0Q0 1 CONTINUE IF (M.NE.0) GO TO 20 D1=1Q0 D2=X DO 5 N=1,NMAX QN=QFLOAT(N) QN1=QFLOAT(N+1) QN2=QFLOAT(2*N+1) D3=(QN2*X*D2-QN*D1)/QN1 DER=QS1*(QN1*QN/QN2)*(-D1+D3) DV1(N)=D2 DV2(N)=DER D1=D2 D2=D3 5 CONTINUE RETURN 20 QMM=QFLOAT(M*M) DO 25 I=1,M I2=I*2 A=A*QSQRT(QFLOAT(I2-1)/QFLOAT(I2))*QS 25 CONTINUE D1=0Q0 D2=A DO 30 N=M,NMAX QN=QFLOAT(N) QN2=QFLOAT(2*N+1) QN1=QFLOAT(N+1) QNM=QSQRT(QN*QN-QMM) QNM1=QSQRT(QN1*QN1-QMM) D3=(QN2*X*D2-QNM*D1)/QNM1 DER=QS1*(-QN1*QNM*D1+QN*QNM1*D3)/QN2 DV1(N)=D2 DV2(N)=DER D1=D2 D2=D3 30 CONTINUE RETURN END C********************************************************************** C * C Calculation of the matrix T = - RG(Q) * (Q**(-1)) * C * C Input infortmation is in COMMON /CTT/ * C Output information is in COMMON /CT/ * C * C********************************************************************** SUBROUTINE TT(NMAX,NCHECK) INCLUDE 'tmq.par.f' IMPLICIT double precision (A-H,O-Z) REAL*16 F(NPN2,NPN2),B(NPN2),WORK(NPN2),COND, * QR(NPN2,NPN2),QI(NPN2,NPN2), * RGQR(NPN2,NPN2),RGQI(NPN2,NPN2), * A(NPN2,NPN2),C(NPN2,NPN2),D(NPN2,NPN2),E(NPN2,NPN2) double precision TR1(NPN2,NPN2),TI1(NPN2,NPN2) COMPLEX*32 ZQ(NPN2,NPN2),ZW(NPN2) COMPLEX*32 ZQR(NPN2,NPN2),ZAFAC(NPN2,NPN2),ZT(NPN2,NPN2), & ZTHETA(NPN2,NPN2) INTEGER IPIV(NPN2),IPVT(NPN2) COMMON /CHOICE/ ICHOICE COMMON /CT/ TR1,TI1 COMMON /CTT/ QR,QI,RGQR,RGQI NDIM=NPN2 NNMAX=2*NMAX IF (ICHOICE.EQ.2) GO TO 5 C Inversion from NAG-LIB or Waterman's method DO I=1,NNMAX DO J=1,NNMAX ZQ(I,J)=QCMPLX(QR(I,J),QI(I,J)) ZAFAC(I,J)=ZQ(I,J) ENDDO ENDDO INFO=0 c CALL F07ARF(NNMAX,NNMAX,ZQ,NPN2,IPIV,INFO) IF (INFO.NE.0) WRITE (6,1100) INFO c CALL F07AWF(NNMAX,ZQ,NPN2,IPIV,ZW,NPN2,INFO) IF (INFO.NE.0) WRITE (6,1100) INFO 1100 FORMAT ('WARNING: info=', I2) DO I=1,NNMAX DO J=1,NNMAX TR=0D0 TI=0D0 DO K=1,NNMAX ARR=RGQR(I,K) ARI=RGQI(I,K) AR=ZQ(K,J) AI=QIMAG(ZQ(K,J)) TR=TR-ARR*AR+ARI*AI TI=TI-ARR*AI-ARI*AR ENDDO TR1(I,J)=TR TI1(I,J)=TI ENDDO ENDDO GOTO 70 C Gaussian elimination 5 DO 10 N1=1,NNMAX DO 10 N2=1,NNMAX F(N1,N2)=QI(N1,N2) 10 CONTINUE IF (NCHECK.EQ.1) THEN CALL INV1(NMAX,F,A) ELSE CALL INVERT(NDIM,NNMAX,F,A,COND,IPVT,WORK,B) ENDIF CALL PROD(QR,A,C,NDIM,NNMAX) CALL PROD(C,QR,D,NDIM,NNMAX) DO 20 N1=1,NNMAX DO 20 N2=1,NNMAX C(N1,N2)=D(N1,N2)+QI(N1,N2) 20 CONTINUE IF (NCHECK.EQ.1) THEN CALL INV1(NMAX,C,QI) ELSE CALL INVERT(NDIM,NNMAX,C,QI,COND,IPVT,WORK,B) ENDIF CALL PROD(A,QR,D,NDIM,NNMAX) CALL PROD(D,QI,QR,NDIM,NNMAX) CALL PROD(RGQR,QR,A,NDIM,NNMAX) CALL PROD(RGQI,QI,C,NDIM,NNMAX) CALL PROD(RGQR,QI,D,NDIM,NNMAX) CALL PROD(RGQI,QR,E,NDIM,NNMAX) DO 30 N1=1,NNMAX DO 30 N2=1,NNMAX TR1(N1,N2)=-A(N1,N2)-C(N1,N2) TI1(N1,N2)= D(N1,N2)-E(N1,N2) 30 CONTINUE 70 RETURN END C********************************************************************** C * C Calculation of the matrix C = A * B . * C All matrices are (N-by-N) * C Declared line dimension of the arrays A,B, and C in the calling * C program is NDIM * C * C********************************************************************** SUBROUTINE PROD(A,B,C,NDIM,N) REAL*16 A(NDIM,N),B(NDIM,N),C(NDIM,N),CIJ DO 10 I=1,N DO 10 J=1,N CIJ=0Q0 DO 5 K=1,N CIJ=CIJ+A(I,K)*B(K,J) 5 CONTINUE C(I,J)=CIJ 10 CONTINUE RETURN END C********************************************************************** SUBROUTINE INV1 (NMAX,F,A) IMPLICIT REAL*16 (A-H,O-Z) INCLUDE 'tmq.par.f' REAL*16 A(NPN2,NPN2),F(NPN2,NPN2),B(NPN1), * WORK(NPN1),Q1(NPN1,NPN1),Q2(NPN1,NPN1), & P1(NPN1,NPN1),P2(NPN1,NPN1) INTEGER IPVT(NPN1),IND1(NPN1),IND2(NPN1) NDIM=NPN1 NN1=(QFLOAT(NMAX)-0.1Q0)*0.5Q0+1Q0 NN2=NMAX-NN1 DO 5 I=1,NMAX IND1(I)=2*I-1 IF(I.GT.NN1) IND1(I)=NMAX+2*(I-NN1) IND2(I)=2*I IF(I.GT.NN2) IND2(I)=NMAX+2*(I-NN2)-1 5 CONTINUE NNMAX=2*NMAX DO 15 I=1,NMAX I1=IND1(I) I2=IND2(I) DO 15 J=1,NMAX J1=IND1(J) J2=IND2(J) Q1(J,I)=F(J1,I1) Q2(J,I)=F(J2,I2) 15 CONTINUE CALL INVERT(NDIM,NMAX,Q1,P1,COND,IPVT,WORK,B) CALL INVERT(NDIM,NMAX,Q2,P2,COND,IPVT,WORK,B) DO 30 I=1,NNMAX DO 30 J=1,NNMAX A(J,I)=0Q0 30 CONTINUE DO 40 I=1,NMAX I1=IND1(I) I2=IND2(I) DO 40 J=1,NMAX J1=IND1(J) J2=IND2(J) A(J1,I1)=P1(J,I) A(J2,I2)=P2(J,I) 40 CONTINUE RETURN END C********************************************************************* C * C Inversion of a square matrix * C * C Input parameters: * C * C A - square (N-by-N) matrix * C NDIM - declared line dimension of the matrix A in the calling * C program * C * C Output information: * C * C X - square (N-by-N) matrix - result of inverting matrix A * C COND - estimate of ill-conditioning of the matrix A * C * C Temporary arrays: IPVT,WORK,B * C * C********************************************************************* SUBROUTINE INVERT (NDIM,N,A,X,COND,IPVT,WORK,B) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 A(NDIM,N),X(NDIM,N),WORK(N),B(N) INTEGER IPVT(N) CALL DECOMP (NDIM,N,A,COND,IPVT,WORK) IF (COND+1Q0.EQ.COND) PRINT 5,COND C IF (COND+1Q0.EQ.COND) STOP 5 FORMAT(' THE MATRIX IS SINGULAR FOR THE GIVEN NUMERICAL ACCURACY ' * ,'COND = ',D12.6) DO 30 I=1,N DO 10 J=1,N B(J)=0Q0 IF (J.EQ.I) B(J)=1Q0 10 CONTINUE CALL SOLVE (NDIM,N,A,B,IPVT) DO 30 J=1,N X(J,I)=B(J) 30 CONTINUE RETURN END C******************************************************************** SUBROUTINE DECOMP (NDIM,N,A,COND,IPVT,WORK) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 A(NDIM,N),COND,WORK(N) INTEGER IPVT(N) IPVT(N)=1 IF(N.EQ.1) GO TO 80 NM1=N-1 ANORM=0Q0 DO 10 J=1,N T=0Q0 DO 5 I=1,N T=T+QABS(A(I,J)) 5 CONTINUE IF (T.GT.ANORM) ANORM=T 10 CONTINUE DO 35 K=1,NM1 KP1=K+1 M=K DO 15 I=KP1,N IF (QABS(A(I,K)).GT.QABS(A(M,K))) M=I 15 CONTINUE IPVT(K)=M IF (M.NE.K) IPVT(N)=-IPVT(N) T=A(M,K) A(M,K)=A(K,K) A(K,K)=T IF (T.EQ.0Q0) GO TO 35 DO 20 I=KP1,N A(I,K)=-A(I,K)/T 20 CONTINUE DO 30 J=KP1,N T=A(M,J) A(M,J)=A(K,J) A(K,J)=T IF (T.EQ.0Q0) GO TO 30 DO 25 I=KP1,N A(I,J)=A(I,J)+A(I,K)*T 25 CONTINUE 30 CONTINUE 35 CONTINUE DO 50 K=1,N T=0Q0 IF (K.EQ.1) GO TO 45 KM1=K-1 DO 40 I=1,KM1 T=T+A(I,K)*WORK(I) 40 CONTINUE 45 EK=1Q0 IF (T.LT.0Q0) EK=-1Q0 IF (A(K,K).EQ.0Q0) GO TO 90 WORK(K)=-(EK+T)/A(K,K) 50 CONTINUE DO 60 KB=1,NM1 K=N-KB T=0Q0 KP1=K+1 DO 55 I=KP1,N T=T+A(I,K)*WORK(K) 55 CONTINUE WORK(K)=T M=IPVT(K) IF (M.EQ.K) GO TO 60 T=WORK(M) WORK(M)=WORK(K) WORK(K)=T 60 CONTINUE YNORM=0Q0 DO 65 I=1,N YNORM=YNORM+QABS(WORK(I)) 65 CONTINUE CALL SOLVE (NDIM,N,A,WORK,IPVT) ZNORM=0Q0 DO 70 I=1,N ZNORM=ZNORM+QABS(WORK(I)) 70 CONTINUE COND=ANORM*ZNORM/YNORM IF (COND.LT.1Q0) COND=1Q0 RETURN 80 COND=1Q0 IF (A(1,1).NE.0Q0) RETURN 90 COND=1Q52 RETURN END C********************************************************************** SUBROUTINE SOLVE (NDIM,N,A,B,IPVT) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 A(NDIM,N),B(N) INTEGER IPVT(N) IF (N.EQ.1) GO TO 50 NM1=N-1 DO 20 K=1,NM1 KP1=K+1 M=IPVT(K) T=B(M) B(M)=B(K) B(K)=T DO 10 I=KP1,N B(I)=B(I)+A(I,K)*T 10 CONTINUE 20 CONTINUE DO 40 KB=1,NM1 KM1=N-KB K=KM1+1 B(K)=B(K)/A(K,K) T=-B(K) DO 30 I=1,KM1 B(I)=B(I)+A(I,K)*T 30 CONTINUE 40 CONTINUE 50 B(1)=B(1)/A(1,1) RETURN END C******************************************************************** C * C Calculation of the expansion coefficients for the (I,Q,U,V) * C representation. * C * C Input parameters: * C * C LAM - wavelength of light * C CSCA - scattering cross section * C TR and TI - elements of the t-matrix. Transferred through * C COMMON /CTM/ * C NMAX - dimension of T(M) matrices * C * C Output infortmation: * C * C ALF1,...,ALF4,BET1,BET2 - expansion coefficients * C LMAX - number of coefficients minus 1 * C * C******************************************************************** SUBROUTINE GSP(NMAX,CSCA,LAM,ALF1,ALF2,ALF3,ALF4,BET1,BET2,LMAX) INCLUDE 'tmq.par.f' IMPLICIT double precision (A-B,D-H,O-Z),COMPLEX*16 (C) REAL*16 LAM double precision SSIGN(900) double precision CSCA,SSI(NPL),SSJ(NPN1), & ALF1(NPL),ALF2(NPL),ALF3(NPL), & ALF4(NPL),BET1(NPL),BET2(NPL), & TR1(NPL1,NPN4),TR2(NPL1,NPN4), & TI1(NPL1,NPN4),TI2(NPL1,NPN4), & G1(NPL1,NPN6),G2(NPL1,NPN6), & AR1(NPN4),AR2(NPN4),AI1(NPN4),AI2(NPN4), & FR(NPN4,NPN4),FI(NPN4,NPN4),FF(NPN4,NPN4) double precision & B1R(NPL1,NPL1,NPN4),B1I(NPL1,NPL1,NPN4), & B2R(NPL1,NPL1,NPN4),B2I(NPL1,NPL1,NPN4), & D1(NPL1,NPN4,NPN4),D2(NPL1,NPN4,NPN4), & D3(NPL1,NPN4,NPN4),D4(NPL1,NPN4,NPN4), & D5R(NPL1,NPN4,NPN4),D5I(NPL1,NPN4,NPN4), & PLUS1(NPN6*NPN4*NPN4*8) double precision & TR11(NPN6,NPN4,NPN4),TR12(NPN6,NPN4,NPN4), & TR21(NPN6,NPN4,NPN4),TR22(NPN6,NPN4,NPN4), & TI11(NPN6,NPN4,NPN4),TI12(NPN6,NPN4,NPN4), & TI21(NPN6,NPN4,NPN4),TI22(NPN6,NPN4,NPN4) COMPLEX*16 CIM(NPN1) COMMON /TMAT/ TR11,TR12,TR21,TR22,TI11,TI12,TI21,TI22 COMMON /CBESS/ B1R,B1I,B2R,B2I COMMON /SS/ SSIGN EQUIVALENCE ( PLUS1(1),TR11(1,1,1) ) EQUIVALENCE (D1(1,1,1),PLUS1(1)), & (D2(1,1,1),PLUS1(NPL1*NPN4*NPN4+1)), & (D3(1,1,1),PLUS1(NPL1*NPN4*NPN4*2+1)), & (D4(1,1,1),PLUS1(NPL1*NPN4*NPN4*3+1)), & (D5R(1,1,1),PLUS1(NPL1*NPN4*NPN4*4+1)) CALL FACT CALL SIGNUM LMAX=2*NMAX L1MAX=LMAX+1 CI=(0D0,1D0) CIM(1)=CI DO 2 I=2,NMAX CIM(I)=CIM(I-1)*CI 2 CONTINUE SSI(1)=1D0 DO 3 I=1,LMAX I1=I+1 SI=DFLOAT(2*I+1) SSI(I1)=SI IF(I.LE.NMAX) SSJ(I)=DSQRT(SI) 3 CONTINUE CI=-CI DO 5 I=1,NMAX SI=SSJ(I) CCI=CIM(I) DO 4 J=1,NMAX SJ=1D0/SSJ(J) CCJ=CIM(J)*SJ/CCI FR(J,I)=CCJ FI(J,I)=CCJ*CI FF(J,I)=SI*SJ 4 CONTINUE 5 CONTINUE NMAX1=NMAX+1 C ***** CALCULATION OF THE ARRAYS B1 AND B2 ***** K1=1 K2=0 K3=0 K4=1 K5=1 K6=2 C PRINT 3300, B1,B2 3300 FORMAT (' b1 and b2') DO 100 N=1,NMAX C ***** CALCULATION OF THE ARRAYS T1 AND T2 ***** DO 10 NN=1,NMAX M1MAX=MIN0(N,NN)+1 DO 6 M1=1,M1MAX M=M1-1 L1=NPN6+M TT1=TR11(M1,N,NN) TT2=TR12(M1,N,NN) TT3=TR21(M1,N,NN) TT4=TR22(M1,N,NN) TT5=TI11(M1,N,NN) TT6=TI12(M1,N,NN) TT7=TI21(M1,N,NN) TT8=TI22(M1,N,NN) T1=TT1+TT2 T2=TT3+TT4 T3=TT5+TT6 T4=TT7+TT8 TR1(L1,NN)=T1+T2 TR2(L1,NN)=T1-T2 TI1(L1,NN)=T3+T4 TI2(L1,NN)=T3-T4 IF(M.EQ.0) GO TO 6 L1=NPN6-M T1=TT1-TT2 T2=TT3-TT4 T3=TT5-TT6 T4=TT7-TT8 TR1(L1,NN)=T1-T2 TR2(L1,NN)=T1+T2 TI1(L1,NN)=T3-T4 TI2(L1,NN)=T3+T4 6 CONTINUE 10 CONTINUE C ***** END OF THE CALCULATION OF THE ARRAYS T1 AND T2 ***** NN1MAX=NMAX1+N DO 40 NN1=1,NN1MAX N1=NN1-1 C ***** CALCULATION OF THE ARRAYS A1 AND A2 ***** CALL CCG(N,N1,NMAX,K1,K2,G1) NNMAX=MIN0(NMAX,N1+N) NNMIN=MAX0(1,IABS(N-N1)) KN=N+NN1 DO 15 NN=NNMIN,NNMAX NNN=NN+1 SIG=SSIGN(KN+NN) M1MAX=MIN0(N,NN)+NPN6 AAR1=0D0 AAR2=0D0 AAI1=0D0 AAI2=0D0 DO 13 M1=NPN6,M1MAX M=M1-NPN6 SSS=G1(M1,NNN) RR1=TR1(M1,NN) RI1=TI1(M1,NN) RR2=TR2(M1,NN) RI2=TI2(M1,NN) IF(M.EQ.0) GO TO 12 M2=NPN6-M RR1=RR1+TR1(M2,NN)*SIG RI1=RI1+TI1(M2,NN)*SIG RR2=RR2+TR2(M2,NN)*SIG RI2=RI2+TI2(M2,NN)*SIG 12 AAR1=AAR1+SSS*RR1 AAI1=AAI1+SSS*RI1 AAR2=AAR2+SSS*RR2 AAI2=AAI2+SSS*RI2 13 CONTINUE XR=FR(NN,N) XI=FI(NN,N) AR1(NN)=AAR1*XR-AAI1*XI AI1(NN)=AAR1*XI+AAI1*XR AR2(NN)=AAR2*XR-AAI2*XI AI2(NN)=AAR2*XI+AAI2*XR 15 CONTINUE C ***** END OF THE CALCULATION OF THE ARRAYS A1 AND A2 **** CALL CCG(N,N1,NMAX,K3,K4,G2) M1=MAX0(-N1+1,-N) M2=MIN0(N1+1,N) M1MAX=M2+NPN6 M1MIN=M1+NPN6 DO 30 M1=M1MIN,M1MAX BBR1=0D0 BBI1=0D0 BBR2=0D0 BBI2=0D0 DO 25 NN=NNMIN,NNMAX NNN=NN+1 SSS=G2(M1,NNN) BBR1=BBR1+SSS*AR1(NN) BBI1=BBI1+SSS*AI1(NN) BBR2=BBR2+SSS*AR2(NN) BBI2=BBI2+SSS*AI2(NN) 25 CONTINUE B1R(NN1,M1,N)=BBR1 B1I(NN1,M1,N)=BBI1 B2R(NN1,M1,N)=BBR2 B2I(NN1,M1,N)=BBI2 30 CONTINUE 40 CONTINUE 100 CONTINUE C ***** END OF THE CALCULATION OF THE ARRAYS B1 AND B2 **** C ***** CALCULATION OF THE ARRAYS D1,D2,D3,D4, AND D5 ***** c PRINT 3301 3301 FORMAT(' d1, d2, ...') DO 200 N=1,NMAX DO 190 NN=1,NMAX M1=MIN0(N,NN) M1MAX=NPN6+M1 M1MIN=NPN6-M1 NN1MAX=NMAX1+MIN0(N,NN) DO 180 M1=M1MIN,M1MAX M=M1-NPN6 NN1MIN=IABS(M-1)+1 DD1=0D0 DD2=0D0 DO 150 NN1=NN1MIN,NN1MAX XX=SSI(NN1) X1=B1R(NN1,M1,N) X2=B1I(NN1,M1,N) X3=B1R(NN1,M1,NN) X4=B1I(NN1,M1,NN) X5=B2R(NN1,M1,N) X6=B2I(NN1,M1,N) X7=B2R(NN1,M1,NN) X8=B2I(NN1,M1,NN) DD1=DD1+XX*(X1*X3+X2*X4) DD2=DD2+XX*(X5*X7+X6*X8) 150 CONTINUE D1(M1,NN,N)=DD1 D2(M1,NN,N)=DD2 180 CONTINUE MMAX=MIN0(N,NN+2) MMIN=MAX0(-N,-NN+2) M1MAX=NPN6+MMAX M1MIN=NPN6+MMIN DO 186 M1=M1MIN,M1MAX M=M1-NPN6 NN1MIN=IABS(M-1)+1 DD3=0D0 DD4=0D0 DD5R=0D0 DD5I=0D0 M2=-M+2+NPN6 DO 183 NN1=NN1MIN,NN1MAX XX=SSI(NN1) X1=B1R(NN1,M1,N) X2=B1I(NN1,M1,N) X3=B2R(NN1,M1,N) X4=B2I(NN1,M1,N) X5=B1R(NN1,M2,NN) X6=B1I(NN1,M2,NN) X7=B2R(NN1,M2,NN) X8=B2I(NN1,M2,NN) DD3=DD3+XX*(X1*X5+X2*X6) DD4=DD4+XX*(X3*X7+X4*X8) DD5R=DD5R+XX*(X3*X5+X4*X6) DD5I=DD5I+XX*(X4*X5-X3*X6) 183 CONTINUE D3(M1,NN,N)=DD3 D4(M1,NN,N)=DD4 D5R(M1,NN,N)=DD5R D5I(M1,NN,N)=DD5I 186 CONTINUE 190 CONTINUE 200 CONTINUE C ***** END OF THE CALCULATION OF THE D-ARRAYS ***** C ***** CALCULATION OF THE EXPANSION COEFFICIENTS ***** C PRINT 3303 3303 FORMAT (' g1, g2, ...') DK=LAM*LAM/(4D0*CSCA*DACOS(-1D0)) DO 300 L1=1,L1MAX G1L=0D0 G2L=0D0 G3L=0D0 G4L=0D0 G5LR=0D0 G5LI=0D0 L=L1-1 SL=SSI(L1)*DK DO 290 N=1,NMAX NNMIN=MAX0(1,IABS(N-L)) NNMAX=MIN0(NMAX,N+L) IF(NNMAX.LT.NNMIN) GO TO 290 CALL CCG(N,L,NMAX,K1,K2,G1) IF(L.GE.2) CALL CCG(N,L,NMAX,K5,K6,G2) NL=N+L DO 280 NN=NNMIN,NNMAX NNN=NN+1 MMAX=MIN0(N,NN) M1MIN=NPN6-MMAX M1MAX=NPN6+MMAX SI=SSIGN(NL+NNN) DM1=0D0 DM2=0D0 DO 270 M1=M1MIN,M1MAX M=M1-NPN6 IF(M.GE.0) SSS1=G1(M1,NNN) IF(M.LT.0) SSS1=G1(NPN6-M,NNN)*SI DM1=DM1+SSS1*D1(M1,NN,N) DM2=DM2+SSS1*D2(M1,NN,N) 270 CONTINUE FFN=FF(NN,N) SSS=G1(NPN6+1,NNN)*FFN G1L=G1L+SSS*DM1 G2L=G2L+SSS*DM2*SI IF(L.LT.2) GO TO 280 DM3=0D0 DM4=0D0 DM5R=0D0 DM5I=0D0 MMAX=MIN0(N,NN+2) MMIN=MAX0(-N,-NN+2) M1MAX=NPN6+MMAX M1MIN=NPN6+MMIN DO 275 M1=M1MIN,M1MAX M=M1-NPN6 SSS1=G2(NPN6-M,NNN) DM3=DM3+SSS1*D3(M1,NN,N) DM4=DM4+SSS1*D4(M1,NN,N) DM5R=DM5R+SSS1*D5R(M1,NN,N) DM5I=DM5I+SSS1*D5I(M1,NN,N) 275 CONTINUE G5LR=G5LR-SSS*DM5R G5LI=G5LI-SSS*DM5I SSS=G2(NPN4,NNN)*FFN G3L=G3L+SSS*DM3 G4L=G4L+SSS*DM4*SI 280 CONTINUE 290 CONTINUE G1L=G1L*SL G2L=G2L*SL G3L=G3L*SL G4L=G4L*SL G5LR=G5LR*SL G5LI=G5LI*SL ALF1(L1)=G1L+G2L ALF2(L1)=G3L+G4L ALF3(L1)=G3L-G4L ALF4(L1)=G1L-G2L BET1(L1)=G5LR*2D0 BET2(L1)=G5LI*2D0 LMAX=L IF(DABS(G1L).LT.1D-6) GO TO 500 300 CONTINUE 500 CONTINUE RETURN END C**************************************************************** C Calculation of the quantities F(N+1)=0.5*ln(n!) C 0.LE.N.LE.899 SUBROUTINE FACT double precision F(900) COMMON /FAC/ F F(1)=0D0 F(2)=0D0 DO 2 I=3,900 I1=I-1 F(I)=F(I1)+0.5D0*DLOG(DFLOAT(I1)) 2 CONTINUE RETURN END C************************************************************ C Calculation of the array SSIGN(N+1)=sign(n) C 0.LE.N.LE.899 SUBROUTINE SIGNUM double precision SSIGN(900) COMMON /SS/ SSIGN SSIGN(1)=1D0 DO 2 N=2,899 SSIGN(N)=-SSIGN(N-1) 2 CONTINUE RETURN END C****************************************************************** C C Calculation of Clebsch-Gordan coefficients C (n,m:n1,m1/nn,mm) C for given n and n1. m1=mm-m, index mm is found from m as C mm=m*k1+k2 C C Input parameters : N,N1,NMAX,K1,K2 C N.LE.NMAX C N.GE.1 C N1.GE.0 C N1.LE.N+NMAX C Output parameters : GG(M+NPN6,NN+1) - array of the corresponding C coefficients C /M/.LE.N C /M1/=/M*(K1-1)+K2/.LE.N1 C NN.LE.MIN(N+N1,NMAX) C NN.GE.MAX(/MM/,/N-N1/) C If K1=1 and K2=0, then 0.LE.M.LE.N SUBROUTINE CCG(N,N1,NMAX,K1,K2,GG) INCLUDE 'tmq.par.f' IMPLICIT double precision (A-H,O-Z) double precision GG(NPL1,NPN6),CD(0:NPN5),CU(0:NPN5) IF(NMAX.LE.NPN4. & AND.0.LE.N1. & AND.N1.LE.NMAX+N. & AND.N.GE.1. & AND.N.LE.NMAX) GO TO 1 PRINT 5001 STOP 5001 FORMAT(' ERROR IN SUBROUTINE CCG') 1 NNF=MIN0(N+N1,NMAX) MIN=NPN6-N MF=NPN6+N IF(K1.EQ.1.AND.K2.EQ.0) MIN=NPN6 DO 100 MIND=MIN,MF M=MIND-NPN6 MM=M*K1+K2 M1=MM-M IF(IABS(M1).GT.N1) GO TO 90 NNL=MAX0(IABS(MM),IABS(N-N1)) IF(NNL.GT.NNF) GO TO 90 NNU=N+N1 NNM=(NNU+NNL)*0.5D0 IF (NNL.EQ.NNU) NNM=NNL CALL CCGIN(N,N1,M,MM,C) CU(NNL)=C IF (NNL.EQ.NNF) GO TO 50 C2=0D0 C1=C DO 7 NN=NNL+1,MIN0(NNM,NNF) A=DFLOAT((NN+MM)*(NN-MM)*(N1-N+NN)) A=A*DFLOAT((N-N1+NN)*(N+N1-NN+1)*(N+N1+NN+1)) A=DFLOAT(4*NN*NN)/A A=A*DFLOAT((2*NN+1)*(2*NN-1)) A=DSQRT(A) B=0.5D0*DFLOAT(M-M1) D=0D0 IF(NN.EQ.1) GO TO 5 B=DFLOAT(2*NN*(NN-1)) B=DFLOAT((2*M-MM)*NN*(NN-1)-MM*N*(N+1)+ & MM*N1*(N1+1))/B D=DFLOAT(4*(NN-1)*(NN-1)) D=D*DFLOAT((2*NN-3)*(2*NN-1)) D=DFLOAT((NN-MM-1)*(NN+MM-1)*(N1-N+NN-1))/D D=D*DFLOAT((N-N1+NN-1)*(N+N1-NN+2)*(N+N1+NN)) D=DSQRT(D) 5 C=A*(B*C1-D*C2) C2=C1 C1=C CU(NN)=C 7 CONTINUE IF (NNF.LE.NNM) GO TO 50 CALL DIRECT(N,M,N1,M1,NNU,MM,C) CD(NNU)=C IF (NNU.EQ.NNM+1) GO TO 50 C2=0D0 C1=C DO 12 NN=NNU-1,NNM+1,-1 A=DFLOAT((NN-MM+1)*(NN+MM+1)*(N1-N+NN+1)) A=A*DFLOAT((N-N1+NN+1)*(N+N1-NN)*(N+N1+NN+2)) A=DFLOAT(4*(NN+1)*(NN+1))/A A=A*DFLOAT((2*NN+1)*(2*NN+3)) A=DSQRT(A) B=DFLOAT(2*(NN+2)*(NN+1)) B=DFLOAT((2*M-MM)*(NN+2)*(NN+1)-MM*N*(N+1) & +MM*N1*(N1+1))/B D=DFLOAT(4*(NN+2)*(NN+2)) D=D*DFLOAT((2*NN+5)*(2*NN+3)) D=DFLOAT((NN+MM+2)*(NN-MM+2)*(N1-N+NN+2))/D D=D*DFLOAT((N-N1+NN+2)*(N+N1-NN-1)*(N+N1+NN+3)) D=DSQRT(D) C=A*(B*C1-D*C2) C2=C1 C1=C CD(NN)=C 12 CONTINUE 50 DO 9 NN=NNL,NNF IF (NN.LE.NNM) GG(MIND,NN+1)=CU(NN) IF (NN.GT.NNM) GG(MIND,NN+1)=CD(NN) c WRITE (6,*) N,M,N1,M1,NN,MM,GG(MIND,NN+1) 9 CONTINUE 90 CONTINUE 100 CONTINUE RETURN END C********************************************************************* SUBROUTINE DIRECT (N,M,N1,M1,NN,MM,C) IMPLICIT double precision (A-H,O-Z) double precision F(900) COMMON /FAC/ F C=F(2*N+1)+F(2*N1+1)+F(N+N1+M+M1+1)+F(N+N1-M-M1+1) C=C-F(2*(N+N1)+1)-F(N+M+1)-F(N-M+1)-F(N1+M1+1)-F(N1-M1+1) C=DEXP(C) RETURN END C********************************************************************* C C Calculation of the Clebcsh-Gordan coefficients C G=(n,m:n1,mm-m/nn,mm) C for given n,n1,m,mm, where NN=MAX(/MM/,/N-N1/) C /M/.LE.N C /MM-M/.LE.N1 C /MM/.LE.N+N1 SUBROUTINE CCGIN(N,N1,M,MM,G) IMPLICIT double precision (A-H,O-Z) double precision F(900),SSIGN(900) COMMON /SS/ SSIGN COMMON /FAC/ F M1=MM-M IF(N.GE.IABS(M). & AND.N1.GE.IABS(M1). & AND.IABS(MM).LE.(N+N1)) GO TO 1 PRINT 5001 STOP 5001 FORMAT(' ERROR IN SUBROUTINE CCGIN') 1 IF (IABS(MM).GT.IABS(N-N1)) GO TO 100 L1=N L2=N1 L3=M IF(N1.LE.N) GO TO 50 K=N N=N1 N1=K K=M M=M1 M1=K 50 N2=N*2 M2=M*2 N12=N1*2 M12=M1*2 G=SSIGN(N1+M1+1) & *DEXP(F(N+M+1)+F(N-M+1)+F(N12+1)+F(N2-N12+2)-F(N2+2) & -F(N1+M1+1)-F(N1-M1+1)-F(N-N1+MM+1)-F(N-N1-MM+1)) N=L1 N1=L2 M=L3 RETURN 100 A=1D0 L1=M L2=MM IF(MM.GE.0) GO TO 150 MM=-MM M=-M M1=-M1 A=SSIGN(MM+N+N1+1) 150 G=A*SSIGN(N+M+1) & *DEXP(F(2*MM+2)+F(N+N1-MM+1)+F(N+M+1)+F(N1+M1+1) & -F(N+N1+MM+2)-F(N-N1+MM+1)-F(-N+N1+MM+1)-F(N-M+1) & -F(N1-M1+1)) M=L1 MM=L2 RETURN END C***************************************************************** SUBROUTINE SAREA (D,RAT) IMPLICIT double precision (A-H,O-Z) IF (D.GE.1) GO TO 10 E=DSQRT(1D0-D*D) R=0.5D0*(D**(2D0/3D0) + D**(-1D0/3D0)*DASIN(E)/E) R=DSQRT(R) RAT=1D0/R RETURN 10 E=DSQRT(1D0-1D0/(D*D)) R=0.25D0*(2D0*D**(2D0/3D0) + D**(-4D0/3D0)*DLOG((1D0+E)/(1D0-E)) & /E) R=DSQRT(R) RAT=1D0/R RETURN END c**************************************************************** SUBROUTINE SURFCH (N,E,RAT) IMPLICIT double precision (A-H,O-Z) double precision X(60),W(60) DN=DFLOAT(N) E2=E*E EN=E*DN NG=60 CALL GAUSS (NG,0,0,X,W) S=0D0 V=0D0 DO 10 I=1,NG XI=X(I) DX=DACOS(XI) DXN=DN*DX DS=DSIN(DX) DSN=DSIN(DXN) DCN=DCOS(DXN) A=1D0+E*DCN A2=A*A ENS=EN*DSN S=S+W(I)*A*DSQRT(A2+ENS*ENS) V=V+W(I)*(DS*A+XI*ENS)*DS*A2 10 CONTINUE RS=DSQRT(S*0.5D0) RV=(V*3D0/4D0)**(1D0/3D0) RAT=RV/RS RETURN END C******************************************************************** SUBROUTINE SAREAC (EPS,RAT) IMPLICIT double precision (A-H,O-Z) RAT=(1.5D0/EPS)**(1D0/3D0) RAT=RAT/DSQRT( (EPS+2D0)/(2D0*EPS) ) RETURN END C******************************************************************** C Computation of R1 and R2 for a power law size distribution with C effective radius A and effective variance B SUBROUTINE POWER (A,B,R1,R2) IMPLICIT double precision (A-H,O-Z) EXTERNAL F COMMON AA,BB AA=A BB=B AX=1D-5 BX=A-1D-5 R1=ZEROIN (AX,BX,F,0D0) R2=(1D0+B)*2D0*A-R1 RETURN END C*********************************************************************** DOUBLE PRECISION FUNCTION ZEROIN (AX,BX,F,TOL) IMPLICIT double precision (A-H,O-Z) EPS=1D0 10 EPS=0.5D0*EPS TOL1=1D0+EPS IF (TOL1.GT.1D0) GO TO 10 15 A=AX B=BX FA=F(A) FB=F(B) 20 C=A FC=FA D=B-A E=D 30 IF (DABS(FC).GE.DABS(FB)) GO TO 40 35 A=B B=C C=A FA=FB FB=FC FC=FA 40 TOL1=2D0*EPS*DABS(B)+0.5D0*TOL XM=0.5D0*(C-B) IF (DABS(XM).LE.TOL1) GO TO 90 44 IF (FB.EQ.0D0) GO TO 90 45 IF (DABS(E).LT.TOL1) GO TO 70 46 IF (DABS(FA).LE.DABS(FB)) GO TO 70 47 IF (A.NE.C) GO TO 50 48 S=FB/FA P=2D0*XM*S Q=1D0-S GO TO 60 50 Q=FA/FC R=FB/FC S=FB/FA P=S*(2D0*XM*Q*(Q-R)-(B-A)*(R-1D0)) Q=(Q-1D0)*(R-1D0)*(S-1D0) 60 IF (P.GT.0D0) Q=-Q P=DABS(P) IF ((2D0*P).GE.(3D0*XM*Q-DABS(TOL1*Q))) GO TO 70 64 IF (P.GE.DABS(0.5D0*E*Q)) GO TO 70 65 E=D D=P/Q GO TO 80 70 D=XM E=D 80 A=B FA=FB IF (DABS(D).GT.TOL1) B=B+D IF (DABS(D).LE.TOL1) B=B+DSIGN(TOL1,XM) FB=F(B) IF ((FB*(FC/DABS(FC))).GT.0D0) GO TO 20 85 GO TO 30 90 ZEROIN=B RETURN END C*********************************************************************** DOUBLE PRECISION FUNCTION F(R1) IMPLICIT double precision (A-H,O-Z) COMMON A,B R2=(1D0+B)*2D0*A-R1 F=(R2-R1)/DLOG(R2/R1)-A RETURN END C********************************************************************** C CALCULATION OF POINTS AND WEIGHTS OF GAUSSIAN QUADRATURE * C FORMULA. IF IND1 = 0 - ON INTERVAL (-1,1), IF IND1 = 1 - ON * C INTERVAL (0,1). IF IND2 = 1 RESULTS ARE PRINTED. * C N - NUMBER OF POINTS * C Z - DIVISION POINTS * C W - WEIGHTS * C********************************************************************** SUBROUTINE GAUSS ( N,IND1,IND2,Z,W ) IMPLICIT double precision (A-H,P-Z) double precision Z(N),W(N) DATA A,B,C /1D0,2D0,3D0/ IND=MOD(N,2) K=N/2+IND F=DFLOAT(N) DO 100 I=1,K M=N+1-I IF(I.EQ.1) X=A-B/((F+A)*F) IF(I.EQ.2) X=(Z(N)-A)*4D0+Z(N) IF(I.EQ.3) X=(Z(N-1)-Z(N))*1.6D0+Z(N-1) IF(I.GT.3) X=(Z(M+1)-Z(M+2))*C+Z(M+3) IF(I.EQ.K.AND.IND.EQ.1) X=0D0 NITER=0 CHECK=1D-16 10 PB=1D0 NITER=NITER+1 IF (NITER.LE.100) GO TO 15 CHECK=CHECK*10D0 15 PC=X DJ=A DO 20 J=2,N DJ=DJ+A PA=PB PB=PC 20 PC=X*PB+(X*PB-PA)*(DJ-A)/DJ PA=A/((PB-X*PC)*F) PB=PA*PC*(A-X*X) X=X-PB IF(DABS(PB).GT.CHECK*DABS(X)) GO TO 10 Z(M)=X W(M)=PA*PA*(A-X*X) IF(IND1.EQ.0) W(M)=B*W(M) IF(I.EQ.K.AND.IND.EQ.1) GO TO 100 Z(I)=-Z(M) W(I)=W(M) 100 CONTINUE IF(IND2.NE.1) GO TO 110 PRINT 1100,N 1100 FORMAT(' *** POINTS AND WEIGHTS OF GAUSSIAN QUADRATURE FORMULA', * ' OF ',I4,'-TH ORDER') DO 105 I=1,K ZZ=-Z(I) 105 PRINT 1200,I,ZZ,I,W(I) 1200 FORMAT(' ',4X,'X(',I4,') = ',F17.14,5X,'W(',I4,') = ',F17.14) GO TO 115 110 CONTINUE C PRINT 1300,N 1300 FORMAT(' GAUSSIAN QUADRATURE FORMULA OF ',I4,'-TH ORDER IS USED') 115 CONTINUE IF(IND1.EQ.0) GO TO 140 DO 120 I=1,N 120 Z(I)=(A+Z(I))/B 140 CONTINUE RETURN END C**************************************************************** SUBROUTINE DISTRB (NNK,YY,WY,NDISTR,AA,BB,GAM,R1,R2,REFF, & VEFF,PI) IMPLICIT double precision (A-H,O-Z) double precision YY(NNK),WY(NNK) IF (NDISTR.EQ.2) GO TO 100 IF (NDISTR.EQ.3) GO TO 200 IF (NDISTR.EQ.4) GO TO 300 IF (NDISTR.EQ.5) GO TO 360 PRINT 1001,AA,BB,GAM 1001 FORMAT('MODIFIED GAMMA DISTRIBUTION, alpha=',F6.4,' r_c=', & F6.4,' gamma=',F6.4) A2=AA/GAM DB=1D0/BB DO 50 I=1,NNK X=YY(I) Y=X**AA X=X*DB Y=Y*DEXP(-A2*(X**GAM)) WY(I)=WY(I)*Y 50 CONTINUE GO TO 400 100 PRINT 1002,AA,BB 1002 FORMAT('LOG-NORMAL DISTRIBUTION, r_g=',F8.4, & ' [ln(sigma_g)]**2=',F6.4) DA=1D0/AA DO 150 I=1,NNK X=YY(I) Y=DLOG(X*DA) Y=DEXP(-Y*Y*0.5D0/BB)/X WY(I)=WY(I)*Y 150 CONTINUE GO TO 400 200 PRINT 1003 1003 FORMAT('POWER LAW DISTRIBUTION OF HANSEN & TRAVIS 1974') DO 250 I=1,NNK X=YY(I) WY(I)=WY(I)/(X*X*X) 250 CONTINUE GO TO 400 300 PRINT 1004,AA,BB 1004 FORMAT ('GAMMA DISTRIBUTION, a=',F8.4,' b=',F6.4) B2=(1D0-3D0*BB)/BB DAB=1D0/(AA*BB) DO 350 I=1,NNK X=YY(I) X=(X**B2)*DEXP(-X*DAB) WY(I)=WY(I)*X 350 CONTINUE GO TO 400 360 PRINT 1005,BB 1005 FORMAT ('MODIFIED POWER LAW DISTRIBUTION, alpha=',D10.4) DO 370 I=1,NNK X=YY(I) IF (X.LE.R1) WY(I)=WY(I) IF (X.GT.R1) WY(I)=WY(I)*(X/R1)**BB 370 CONTINUE 400 CONTINUE SUM=0D0 DO 450 I=1,NNK SUM=SUM+WY(I) 450 CONTINUE SUM=1D0/SUM DO 500 I=1,NNK WY(I)=WY(I)*SUM 500 CONTINUE G=0D0 DO 550 I=1,NNK X=YY(I) G=G+X*X*WY(I) 550 CONTINUE REFF=0D0 DO 600 I=1,NNK X=YY(I) REFF=REFF+X*X*X*WY(I) 600 CONTINUE REFF=REFF/G VEFF=0D0 DO 650 I=1,NNK X=YY(I) XI=X-REFF VEFF=VEFF+XI*XI*X*X*WY(I) 650 CONTINUE VEFF=VEFF/(G*REFF*REFF) RETURN END C************************************************************* SUBROUTINE HOVENR(L1,A1,A2,A3,A4,B1,B2) IMPLICIT double precision (A-H,O-Z) double precision + A1(L1),A2(L1),A3(L1),A4(L1),B1(L1),B2(L1) DO 100 L=1,L1 KONTR=1 LL=L-1 DL=DFLOAT(LL)*2D0+1D0 DDL=DL*0.48D0 AA1=A1(L) AA2=A2(L) AA3=A3(L) AA4=A4(L) BB1=B1(L) BB2=B2(L) IF(LL.GE.1.AND.DABS(AA1).GE.DL) KONTR=2 IF(DABS(AA2).GE.DL) KONTR=2 IF(DABS(AA3).GE.DL) KONTR=2 IF(DABS(AA4).GE.DL) KONTR=2 IF(DABS(BB1).GE.DDL) KONTR=2 IF(DABS(BB2).GE.DDL) KONTR=2 IF(KONTR.EQ.2) PRINT 3000,LL C=-0.1D0 DO 50 I=1,11 C=C+0.1D0 CC=C*C C1=CC*BB2*BB2 C2=C*AA4 C3=C*AA3 IF((DL-C*AA1)*(DL-C*AA2)-CC*BB1*BB1.LE.-1D-4) KONTR=2 IF((DL-C2)*(DL-C3)+C1.LE.-1D-4) KONTR=2 IF((DL+C2)*(DL-C3)-C1.LE.-1D-4) KONTR=2 IF((DL-C2)*(DL+C3)-C1.LE.-1D-4) KONTR=2 IF(KONTR.EQ.2) PRINT 4000,LL,C 50 CONTINUE 100 CONTINUE IF(KONTR.EQ.1) PRINT 2000 2000 FORMAT('TEST OF VAN DER MEE & HOVENIER IS SATISFIED') 3000 FORMAT('TEST OF VAN DER MEE & HOVENIER IS NOT SATISFIED, L=',I3) 4000 FORMAT('TEST OF VAN DER MEE & HOVENIER IS NOT SATISFIED, L=',I3, & ' A=',D9.2) RETURN END C**************************************************************** C CALCULATION OF THE SCATTERING MATRIX FOR GIVEN EXPANSION C COEFFICIENTS C A1,...,B2 - EXPANSION COEFFICIENTS C LMAX - NUMBER OF COEFFICIENTS MINUS 1 C N - NUMBER OF SCATTERING ANGLES C THE CORRESPONDING SCATTERING ANGLES ARE GIVEN BY C 180*(I-1)/(N-1) (DEGREES), WHERE I NUMBERS THE ANGLES SUBROUTINE MATR(A1,A2,A3,A4,B1,B2,LMAX,NPNA) INCLUDE 'tmq.par.f' IMPLICIT double precision (A-H,O-Z) double precision + A1(NPL),A2(NPL),A3(NPL),A4(NPL),B1(NPL),B2(NPL) N=NPNA DN=1D0/DFLOAT(N-1) DA=DACOS(-1D0)*DN DB=180D0*DN L1MAX=LMAX+1 PRINT 1000 1000 FORMAT(' ') PRINT 1001 1001 FORMAT(' ',2X,'S',6X,'ALPHA1',6X,'ALPHA2',6X,'ALPHA3', & 6X,'ALPHA4',7X,'BETA1',7X,'BETA2') DO 10 L1=1,L1MAX L=L1-1 PRINT 1002,L,A1(L1),A2(L1),A3(L1),A4(L1),B1(L1),B2(L1) 10 CONTINUE 1002 FORMAT(' ',I3,6F12.5) TB=-DB TAA=-DA PRINT 1000 PRINT 1003 1003 FORMAT(' ',5X,'<',8X,'F11',8X,'F22',8X,'F33', & 8X,'F44',8X,'F12',8X,'F34') D6=DSQRT(6D0)*0.25D0 DO 500 I1=1,N TAA=TAA+DA TB=TB+DB U=DCOS(TAA) F11=0D0 F2=0D0 F3=0D0 F44=0D0 F12=0D0 F34=0D0 P1=0D0 P2=0D0 P3=0D0 P4=0D0 PP1=1D0 PP2=0.25D0*(1D0+U)*(1D0+U) PP3=0.25D0*(1D0-U)*(1D0-U) PP4=D6*(U*U-1D0) DO 400 L1=1,L1MAX L=L1-1 DL=DFLOAT(L) DL1=DFLOAT(L1) F11=F11+A1(L1)*PP1 F44=F44+A4(L1)*PP1 IF(L.EQ.LMAX) GO TO 350 PL1=DFLOAT(2*L+1) P=(PL1*U*PP1-DL*P1)/DL1 P1=PP1 PP1=P 350 IF(L.LT.2) GO TO 400 F2=F2+(A2(L1)+A3(L1))*PP2 F3=F3+(A2(L1)-A3(L1))*PP3 F12=F12+B1(L1)*PP4 F34=F34+B2(L1)*PP4 IF(L.EQ.LMAX) GO TO 400 PL2=DFLOAT(L*L1)*U PL3=DFLOAT(L1*(L*L-4)) PL4=1D0/DFLOAT(L*(L1*L1-4)) P=(PL1*(PL2-4D0)*PP2-PL3*P2)*PL4 P2=PP2 PP2=P P=(PL1*(PL2+4D0)*PP3-PL3*P3)*PL4 P3=PP3 PP3=P P=(PL1*U*PP4-DSQRT(DFLOAT(L*L-4))*P4)/DSQRT(DFLOAT(L1*L1-4)) P4=PP4 PP4=P 400 CONTINUE F22=(F2+F3)*0.5D0 F33=(F2-F3)*0.5D0 C F22=F22/F11 C F33=F33/F11 C F44=F44/F11 C F12=-F12/F11 C F34=F34/F11 PRINT 1004,TB,F11,F22,F33,F44,F12,F34 500 CONTINUE PRINT 1000 1004 FORMAT(' ',F6.2,6F11.4) RETURN END