C
C  Fortran code gmm03sq.f for calculation of elastic radiative scattering 
C  by a single or an ensemble of variously shaped small particles in a 
C  single, fixed orientation
C  In all aspects and usage this gmm03sq.f is exactly the same as the 
C  gmm03s.f code except that the double precision subroutine tm0d.f in 
C  the gmm03s.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 gmm03sq.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 gmm03sq.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 
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      gmm03sq.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      crgmm03sq.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 gmm03sq
      implicit double precision (a-h,o-z)
      include 'gmm01f.par'
      include 'tmq.par.f'
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 for "gmm01f.par":   
C     parameter (nLp=1000,np=10) 
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  **********  Computer memory required by this code is at the level of
C                                np^4 + nLp*np^3
C  This gmm03sq.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 "tmq.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-----------------------------------------------------------------------
      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(2,2,nmp,nmp),
     +   tbar0(nLp,2,2,0:np,np,np),tbars(2,2,0:np,np,np),
     +   ealp(nLp,-2*np:2*np),ealps(-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/tmtrot/ealps,tbars,tbar
      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               "tmq.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'
         write(6,*)
     +      'rotational symmetry of nonspherical component particles'
         write(6,*)
     +      'are randomly picked up by this code  ******'
      endif
      if(idbeta.eq.1) then
         write(6,*)
     +      '******  Note that the zenith angle of the axis of'
         write(6,*)
     +      'rotational symmetry of nonspherical component particles'
         write(6,*)
     +      'are 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<fint, interaction between the pair is 
C  considered to be negligible and no interaction will be calculated 
C  between the two particles. When fint=0, no particle is excluded in 
C  the interaction calculations.
C-----------------------------------------------------------------------
      fint=0.d0
      if(fint.lt.0.d0.or.fint.gt.1.d0) then
         fint=0.02d0
         write(6,'(a37)') 'Interaction index: using default 0.02'
      else
         write(6,'(a18,f7.3)') 'Interaction index:',fint
      endif
      READ(1,*) sang,pang
C-----------------------------------------------------------------------
C  sang -- the scattering angle interval for the output of the amplitude 
C       scattering matrix and Mueller matrix
C       example:   
C       when sang=1, the results in the output files of "amp.forw.out", 
C       "amp.back.out", "mue.forw.out", and "mue.back.out" will be 
C       written for every degree of scattering angle in the range of 
C       [0,90] or [90,180]
C       ******  Note it is better to have 90/sang as an integer  ******
C  pang -- the azimuth angle interval in the output of the 
C       two-dimensional amplitude and Mueller matrices in the range 
C       of [0,360]
C       when pang=0, calculates only the scattering plane of phi=0, 
C       i.e., the x-z plane
C-----------------------------------------------------------------------
      write(6,'(a32,f7.3)') 
     +     'scat.-angle-interval in output: ',sang
      if(sang.le.0.d0) sang=1.d0
      nang=90.d0/sang+1.d0
      nang2=2*nang-1
      if(nang2.gt.nangmax) then
         write(6,*) 'sang too small'
         write(6,*) 
     +      'please increase sang in the input file gmm03s.in'
         write(6,*) '  and try again, or'
         write(6,*) 
     +      'increase nangmax in the parameter line of the'
         write(6,*)  'main code, recompile, then try again'
	 stop
      endif
      write(6,'(a,f8.3)') 
     +   'azimuth-angle-interval in Mueller matrix output: ',pang
      if(pang.lt.0.0001d0) then
         npng=1
      else
         npng=360.0d0/pang
      endif
      close(1)
      write(6,'(/)')
      OPEN(UNIT=2,FILE=FLNAME,STATUS='OLD')
      READ(2,*) w,irat
C-----------------------------------------------------------------------
C  w -- incident wavelength
C  irat -- the same as the parameter RAT in Mishchenko's code
C     when irat=1, particle size is specified in terms of the 
C     equal-volume-sphere radius, otherwise, particle size is specified 
C     in terms of the equal-surface-area-sphere radius
C-----------------------------------------------------------------------
      READ(2,*) nL
C-----------------------------------------------------------------------
C  nL -- number of spheres in the aggregate
C-----------------------------------------------------------------------
      if(nL.gt.nLp) then
         write(6,*) 'Parameter nLp too small, must be >', nL 
         write(6,*) 
     +      'Change nLp in gmm01f.par, recompile, then try again'
         stop
      endif
      if(nL.eq.1) 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
         if(idshp(i).ne.0) then
            r0(7,i)=0.d0
            r0(8,i)=0.d0
            r0(9,i)=0.d0
         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='gmm03sq.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 m=0,np   
            do n=1,np
               do v=1,np
	          do ip=1,2
		     do iq=1,2
                        tbar0(i,ip,iq,m,n,v)=0.d0
                     enddo
		  enddo
	       enddo
            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)
			   do m=0,min(n,v)
                              do ip=1,2
                                 do iq=1,2
			            A=tbar0(j,ip,iq,m,n,v)
			            tbar0(i,ip,iq,m,n,v)=A
                                 enddo
                              enddo
			   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 n=1,nmax(i)
               do m=0,n
                  tbar0(i,1,1,m,n,n)=an(n)
                  tbar0(i,2,2,m,n,n)=bn(n)
	       enddo   
            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 tmq.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)
               do m=0,nsmall
                  B=dcmplx(RT11(m+1,n,v),IT11(m+1,n,v))
                  tbar0(i,1,1,m,n,v)=A*B
                  B=dcmplx(RT12(m+1,n,v),IT12(m+1,n,v))
                  tbar0(i,1,2,m,n,v)=A*B
                  B=dcmplx(RT21(m+1,n,v),IT21(m+1,n,v))
                  tbar0(i,2,1,m,n,v)=A*B
                  B=dcmplx(RT22(m+1,n,v),IT22(m+1,n,v))
                  tbar0(i,2,2,m,n,v)=A*B
               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)
         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 170
         sa=dsin(alph)
         ca=dcos(alph)
         A=dcmplx(ca,sa)
         n1=nmax(i)
	 ealp(i,0)=1.d0
         do m=1,2*n1
            ealp(i,m)=A**m
	    ealp(i,-m)=dconjg(ealp(i,m))
         enddo
 170  enddo
      cextx=0.d0
      cexty=0.d0
      cabsx=0.d0
      cabsy=0.d0
      cscax=0.d0
      cscay=0.d0
      cprx=0.d0
      cpry=0.d0
      cbakx=0.d0
      cbaky=0.d0
      do i=1,nL
         cextxi(i)=0.d0
         cextyi(i)=0.d0
         cabsxi(i)=0.d0
         cabsyi(i)=0.d0
         cscaxi(i)=0.d0
         cscayi(i)=0.d0
         cprxi(i)=0.d0
         cpryi(i)=0.d0
      enddo
      do i=1,nang2
         i11(i)=0.d0
         i21(i)=0.d0
         i22(i)=0.d0
         i12(i)=0.d0
      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  the Gaunt coefficients are required in the evaluation of vector 
C  translation coefficients; formulations used in this code for vector 
C  translation coefficients are from Cruzan, Q. Appl. Math. 20, 33 
C  (1962) and  Xu, J. Comput. 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
         N1=nmax(i)
         if(idshp(i).eq.0) goto 202 
         sa=shp(2,i)         
         sb=shp(3,i)
         if(sa.eq.0.d0.and.sb.eq.0.d0) goto 202
         if(i.eq.1) goto 201
	 ii=idshp(i)
         ij=idshp(i-1)
	 s=shp(1,i)
	 t=shp(1,i-1)
         if(ii.eq.ij.and.s.eq.t) then
            if(x(i).eq.x(i-1).and.ref(i).eq.ref(i-1)) then
               if(xc(i).eq.xc(i-1).and.refc(i).eq.refc(i-1)) then
                  ca=shp(2,i-1)
                  cb=shp(3,i-1)
                  if(sa.eq.ca.and.sb.eq.cb) goto 202
               endif
            endif
         endif         
 201     do m=-2*N1,2*N1
            ealps(m)=ealp(i,m)
	 enddo         
         do m=0,N1
            do n=1,N1
               do v=1,N1
                  do ip=1,2
                     do iq=1,2
                        A=tbar0(i,ip,iq,m,n,v)
                        tbars(ip,iq,m,n,v)=A
                     enddo
                  enddo
               enddo
            enddo
         enddo
         call tmxbar(N1,sa,sb,pih)
 202     cz=dcos(k*r0(3,i))
         sz=dsin(k*r0(3,i))
         cmz=0.5d0*dcmplx(cz,sz)
         do 21 n=1,N1
            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=tbar0(i,1,1,0,n,n)
	       as(i,imn)=B*A
	       as(i,imn-2)=B*p0(i,imn-2)
	       B=tbar0(i,2,2,0,n,n)
	       bs(i,imn)=B*A
	       bs(i,imn-2)=B*q0(i,imn-2)
	    else
               if(sa.eq.0.d0.and.sb.eq.0.d0) then
                  do v=1,N1
	             iuv=v*(v+1)+1
                     B=tbar0(i,1,1,1,v,n)
                     B=B+tbar0(i,1,2,1,v,n)
                     as(i,iuv)=as(i,iuv)+B*A
                     B=tbar0(i,2,1,1,v,n)
                     B=B+tbar0(i,2,2,1,v,n)
                     bs(i,iuv)=bs(i,iuv)+B*A
                     B=tbar0(i,1,1,1,v,n)
                     Aj=B*p0(i,imn-2)
                     B=-tbar0(i,1,2,1,v,n)
                     Aj=Aj+B*q0(i,imn-2)
                     as(i,iuv-2)=as(i,iuv-2)+Aj
                     B=-tbar0(i,2,1,1,v,n)
                     Aj=B*p0(i,imn-2)
                     B=tbar0(i,2,2,1,v,n)
                     Aj=Aj+B*q0(i,imn-2)
                     bs(i,iuv-2)=bs(i,iuv-2)+Aj
                  enddo
               else
                  do iuv=1,uvmax(i)
                     B=tbar(1,1,iuv,imn)
                     B=B+tbar(1,2,iuv,imn)
                     A0=B*A
                     as(i,iuv)=as(i,iuv)+A0
                     B=tbar(1,1,iuv,imn-2)
                     A0=B*p0(i,imn-2)
                     B=tbar(1,2,iuv,imn-2)
                     A0=A0+B*q0(i,imn-2)
                     as(i,iuv)=as(i,iuv)+A0
                     B=tbar(2,1,iuv,imn)
                     B=B+tbar(2,2,iuv,imn)
                     A0=B*A
                     bs(i,iuv)=bs(i,iuv)+A0
                     B=tbar(2,1,iuv,imn-2)
                     A0=B*p0(i,imn-2)
                     B=tbar(2,2,iuv,imn-2)
                     A0=A0+B*q0(i,imn-2)
                     bs(i,iuv)=bs(i,iuv)+A0		  
                  enddo	       
               endif
            endif
 21      continue
         if(idshp(i).eq.0.or.(sa.eq.0.d0.and.sb.eq.0.d0)) then
            do n=1,N1
               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
         N1=nmax(i)
         if(idshp(i).eq.0) goto 602 
         sa=shp(2,i)
         sb=shp(3,i)
         if(sa.eq.0.d0.and.sb.eq.0.d0) goto 602
         if(i.eq.1) goto 601
	 ii=idshp(i)
         ij=idshp(i-1)
	 s=shp(1,i)
	 t=shp(1,i-1)
         if(ii.eq.ij.and.s.eq.t) then
            if(x(i).eq.x(i-1).and.ref(i).eq.ref(i-1)) then
               if(xc(i).eq.xc(i-1).and.refc(i).eq.refc(i-1)) then
                  ca=shp(2,i-1)
                  cb=shp(3,i-1)
                  if(sa.eq.ca.and.sb.eq.cb) goto 602
               endif
            endif
         endif         
 601     do m=-2*N1,2*N1
            ealps(m)=ealp(i,m)
	 enddo         
         do m=0,N1
            do n=1,N1
               do v=1,N1
                  do ip=1,2
                     do iq=1,2
                        A=tbar0(i,ip,iq,m,n,v)
                        tbars(ip,iq,m,n,v)=A
                     enddo
                  enddo
               enddo
            enddo
         enddo
         call tmxbar(N1,sa,sb,pih)
 602     c1i(i)=0.d0
	 if(idshp(i).eq.0) then
            do imn=1,uvmax(i)
               n=dsqrt(dble(imn))
               s2x=tbar0(i,1,1,0,n,n)
               s4x=tbar0(i,2,2,0,n,n)
               as1(i,imn)=-s2x*as00(i,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(sa.eq.0.d0.and.sb.eq.0.d0) then
               do imn=1,uvmax(i)
                  n=dsqrt(dble(imn))
                  m=imn-n*(n+1)
                  m1=iabs(m)
                  s3y=0.d0
                  s1y=0.d0
                  do v=1,nmax(i)
                     if(m1.le.v) then
                        iuv=m+v*(v+1)
                        s2x=tbar0(i,1,1,m1,n,v)
                        s3y=s3y+s2x*as00(i,iuv)
                        s2x=tbar0(i,1,2,m1,n,v)
                        if(m.lt.0) s2x=-s2x
                        s3y=s3y+s2x*bs00(i,iuv)
                        s4x=tbar0(i,2,1,m1,n,v)
                        if(m.lt.0) s4x=-s4x
                        s1y=s1y+s4x*as00(i,iuv)
                        s4x=tbar0(i,2,2,m1,n,v)
                        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(1,1,imn,iuv)
                     s3y=s3y+s2x*as00(i,iuv)
                     s2x=tbar(1,2,imn,iuv)
                     s3y=s3y+s2x*bs00(i,iuv)
                     s4x=tbar(2,1,imn,iuv)
                     s1y=s1y+s4x*as00(i,iuv)
                     s4x=tbar(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
         N1=nmax(i)
         if(idshp(i).eq.0) goto 6212 
         sa=shp(2,i)
         sb=shp(3,i)
         if(sa.eq.0.d0.and.sb.eq.0.d0) goto 6212
         if(i.eq.1) goto 6211
	 ii=idshp(i)
         ij=idshp(i-1)
	 s=shp(1,i)
	 t=shp(1,i-1)
         if(ii.eq.ij.and.s.eq.t) then
            if(x(i).eq.x(i-1).and.ref(i).eq.ref(i-1)) then
               if(xc(i).eq.xc(i-1).and.refc(i).eq.refc(i-1)) then
                  ca=shp(2,i-1)
                  cb=shp(3,i-1)
                  if(sa.eq.ca.and.sb.eq.cb) goto 6212
               endif
            endif
         endif         
6211     do m=-2*N1,2*N1
            ealps(m)=ealp(i,m)
	 enddo         
         do m=0,N1
            do n=1,N1
               do v=1,N1
                  do ip=1,2
                     do iq=1,2
                        A=tbar0(i,ip,iq,m,n,v)
                        tbars(ip,iq,m,n,v)=A
                     enddo
                  enddo
               enddo
            enddo
         enddo
         call tmxbar(N1,sa,sb,pih)
6212     if(idshp(i).eq.0) then
            do imn=1,uvmax(i)
               n=dsqrt(dble(imn))
               s2x=tbar0(i,1,1,0,n,n)
               asv(i,imn)=s2x*as00(i,imn)
               s4x=tbar0(i,2,2,0,n,n)
               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(sa.eq.0.d0.and.sb.eq.0.d0) then
               do imn=1,uvmax(i)
                  n=dsqrt(dble(imn))
                  m=imn-n*(n+1)
                  m1=iabs(m)
                  s3y=0.d0
                  s1y=0.d0
                  do v=1,nmax(i)
                     if(iabs(m).le.v) then
                        iuv=m+v*(v+1)
                        s2x=tbar0(i,1,1,m1,n,v)
                        s3y=s3y+s2x*as00(i,iuv)
                        s2x=tbar0(i,1,2,m1,n,v)
                        if(m.lt.0) s2x=-s2x
                        s3y=s3y+s2x*bs00(i,iuv)
                        s4x=tbar0(i,2,1,m1,n,v)
                        if(m.lt.0) s4x=-s4x
                        s1y=s1y+s4x*as00(i,iuv)
                        s4x=tbar0(i,2,2,m1,n,v)
                        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(1,1,imn,iuv)
                     s3y=s3y+s2x*as00(i,iuv)
                     s2x=tbar(1,2,imn,iuv)
                     s3y=s3y+s2x*bs00(i,iuv)
                     s4x=tbar(2,1,imn,iuv)
                     s1y=s1y+s4x*as00(i,iuv)
                     s4x=tbar(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
         N1=nmax(i)
         if(idshp(i).eq.0) goto 6242 
         sa=shp(2,i)
         sb=shp(3,i)
         if(sa.eq.0.d0.and.sb.eq.0.d0) goto 6242
         if(i.eq.1) goto 6241
	 ii=idshp(i)
         ij=idshp(i-1)
	 s=shp(1,i)
	 t=shp(1,i-1)
         if(ii.eq.ij.and.s.eq.t) then
            if(x(i).eq.x(i-1).and.ref(i).eq.ref(i-1)) then
               if(xc(i).eq.xc(i-1).and.refc(i).eq.refc(i-1)) then
                  ca=shp(2,i-1)
                  cb=shp(3,i-1)
                  if(sa.eq.ca.and.sb.eq.cb) goto 6242
               endif
            endif
         endif         
6241     do m=-2*N1,2*N1
            ealps(m)=ealp(i,m)
	 enddo         
         do m=0,N1
            do n=1,N1
               do v=1,N1
                  do ip=1,2
                     do iq=1,2
                        A=tbar0(i,ip,iq,m,n,v)
                        tbars(ip,iq,m,n,v)=A
                     enddo
                  enddo
               enddo
            enddo
         enddo
         call tmxbar(N1,sa,sb,pih)
6242     if(idshp(i).eq.0) then
            do imn=1,uvmax(i)
               n=dsqrt(dble(imn))
               s2x=tbar0(i,1,1,0,n,n)
               ast(i,imn)=s2x*as00(i,imn)
               s4x=tbar0(i,2,2,0,n,n)
               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(sa.eq.0.d0.and.sb.eq.0.d0) then
               do imn=1,uvmax(i)
                  n=dsqrt(dble(imn))
                  m=imn-n*(n+1)
                  m1=iabs(m)
                  s3y=0.d0
                  s1y=0.d0
                  do v=1,nmax(i)
                     if(iabs(m).le.v) then
                        iuv=m+v*(v+1)
                        s2x=tbar0(i,1,1,m1,n,v)
                        s3y=s3y+s2x*as00(i,iuv)
                        s2x=tbar0(i,1,2,m1,n,v)
                        if(m.lt.0) s2x=-s2x
                        s3y=s3y+s2x*bs00(i,iuv)
                        s4x=tbar0(i,2,1,m1,n,v)
                        if(m.lt.0) s4x=-s4x
                        s1y=s1y+s4x*as00(i,iuv)
                        s4x=tbar0(i,2,2,m1,n,v)
                        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(1,1,imn,iuv)
                     s3y=s3y+s2x*as00(i,iuv)
                     s2x=tbar(1,2,imn,iuv)
                     s3y=s3y+s2x*bs00(i,iuv)
                     s4x=tbar(2,1,imn,iuv)
                     s1y=s1y+s4x*as00(i,iuv)
                     s4x=tbar(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 by
C  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 [see, for
c  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','<cos(theta)>'
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','<cos(theta)>'
      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','<cos(theta)>'
      else
	 write(6,221)
     +    'Qexts','Qabss','Qscas','Qbaks','Qprs','<cos(theta)>'
      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 gmm03sq.f code for'
      write(12,*) '  a single or ensembles of variouly shaped small'
      write(12,*) '  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','<cos(theta)>'
      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','<cos(theta)>'
      else
         write(12,221)
     +    'Qexts','Qabss','Qscas','Qbaks','Qprs','<cos(theta)>'
      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))
            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

      subroutine tmxbar(nmax,shp2,shp3,pih)
      implicit double precision (a-h,o-z)
      include 'gmm01f.par'
      parameter (nmp=np*(np+2))
      integer u,v,uvmax
      complex*16 A,B,Bj,bar(2,2),ealps(-2*np:2*np),
     +   tbars(2,2,0:np,np,np),tbar(2,2,nmp,nmp)
      common/rot/bcof(0:np+2),dc(-np:np,0:nmp)
      common/tmtrot/ealps,tbars,tbar
      uvmax=nmax*(nmax+2)
      do imn=1,uvmax
         do iuv=1,uvmax
            do ip=1,2
               do iq=1,2
                  tbar(ip,iq,imn,iuv)=0.d0
               enddo
            enddo
         enddo
      enddo
      s=shp2*pih/90.d0
      s=dabs(s)
      t=shp3*pih/90.d0
      t0=dabs(t)
      if(s.lt.1.d-10.and.t0.lt.1.d-10) return
      cb=dcos(t)
      call rotcoef(cb,nmax)
      do m=-nmax,nmax
         do u=-nmax,nmax
            A=ealps(u-m)
            do 10 n=1,nmax
               if(iabs(m).gt.n) go to 10
               in0=n*(n+1)
	       imn=in0+m
               do 12 v=1,nmax
                  if(iabs(u).gt.v) go to 12
                  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
                           if(is.lt.0) then
                              temp=(-1)**(ip+iq)
                              Bj=temp*tbars(ip,iq,-is,n,v)
                           else
                              Bj=tbars(ip,iq,is,n,v)
                           endif
                           bar(ip,iq)=bar(ip,iq)+B*Bj
                        enddo
                     enddo	 
                  enddo      
                  do ip=1,2
                     do iq=1,2
                        tbar(ip,iq,imn,iuv)=bar(ip,iq)
                     enddo
                  enddo
 12            continue
 10         continue
         enddo
      enddo
      return
      end

C -------------------------------------------------------------------
C
C  the following is a part of Mishchenko's public domain code 
C  "tmq.new.f" available online at http://www.giss.nasa.gov/~crmin. 
C  For detailed information on the formulation and algorithms, please 
C  see the original code tmq.new.f by Mishchenko.
C  Yu-lin Xu  
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,'<COS>=',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
