C
C  Fortran code gmm01TrA.f for calculating scattering of electromagnetic
C  radiation by a randomly oriented aggregate of homogeneous spheres 
C  (modified from gmm01f.f) 
C  Yu-lin Xu 
C  7/2002: released to public
C
C-----------------------------------------------------------------------
C
C  For questions/comments/suggestions/bugs/problems please contact 
C  Yu-lin Xu at xu@astro.ufl.edu or yxu@erc.ufl.edu
C
C-----------------------------------------------------------------------
C
C  October 2002:  Thanks Piotr Flatau for debugging. 
C
C-----------------------------------------------------------------------
       PROGRAM gmm01TrA
       implicit double precision (a-h,o-z)
       include 'gmm01f.par'
C-----------------------------------------------------------------------
C  Two parameters in "gmm01f.par": nLp, np
C  nLp --- the maximum number of spheres, must be equal to or greater 
C          than the sphere-number in an aggregate actually calculated  
C  np  --- the maximum scattering order in the incident and scattered
C          field expansions, must be equal to or greater than the 
C          highest scattering order required in actual calculations;
C          usually, this can be estimated by Wiscombe's criterion  
C          regarding the largest individual sphere in an aggregate
C  an example of "gmm01f.par":   parameter (nLp=6,np=10)
C  
C
C  ****** This code is more memory-demanding than the single-orientation 
C         code gmm01f.f. The computer memory required is at the level of 
C                         (nLp*np^2)^2.
C   This code is written in double precision arithmetic. An individual  
C   sphere can have a size parameter way beyond ~200. There is no limit  
C   to the overall dimension of an aggregate. The largest individual  
C   sphere size and the maximum number of spheres that an aggregate can  
C   have depends on the availability of computer memory.  ****** 
C
C-----------------------------------------------------------------------
C
C  This code requires two input files: gmmTr.in and a user-specified 
C  datafile for aggregate configuration
C
C  ****** An example for the input and output files ******
C
C  (gmmTr.in)
C       gold2.k
C       0
C       1.d-12  200
C       -1
C       0  90.
C
C  (gold2.k)
C        389.3652205    
C        2
C        0.   0.    0.    6.2   0.52362125   1.63011385
C        0.   0.   14.8   6.2   0.52362125   1.63011385
C
C  (gmm01TrA.out)
C gmm01TrA.out    input file:  gold2.k           xv:   0.126  xs:   0.141
C
C    <Cext>       <Cabs>       <Csca>       <Cbak>       <Cpr>     <cos(theta)>
C  0.17001E+03  0.16940E+03  0.61371E+00  0.87511E+00  0.17001E+03  0.39623E-02
C   <Qextv>      <Qabsv>      <Qscav>      <Qbakv>      <Qprv>     <cos(theta)>
C  0.88687E+00  0.88367E+00  0.32014E-02  0.45650E-02  0.12685E-04  0.39623E-02
C   <Qexts>      <Qabss>      <Qscas>      <Qbaks>      <Qprs>     <cos(theta)>
C  0.70391E+00  0.70137E+00  0.25410E-02  0.36233E-02  0.10068E-04  0.39623E-02
C
C  s.a.    <total>    <pol.>    <S1*S1>      <S4*S4>      <S3*S3>      <S2*S2>
C   0.0  0.37795E-04  0.0000  0.18527E-04  0.37062E-06  0.37062E-06  0.18527E-04
C  90.0  0.19077E-04  0.9412  0.18330E-04  0.37284E-06  0.18654E-06  0.18769E-06
C 180.0  0.37012E-04  0.0000  0.18134E-04  0.37157E-06  0.37157E-06  0.18134E-04
C
C Scattering matrix (4X4 for each scattering angle):
C    0.0   0.1889741E-04   0.0000000E+00   0.3855997E-20   0.1214528E-19
C          0.0000000E+00   0.1815616E-04   0.3855997E-20   0.1214528E-19
C         -0.9639992E-21  -0.2494377E-19   0.1815616E-04  -0.1623600E-21
C          0.3036320E-20   0.6408125E-22   0.4531011E-21   0.1889741E-04
C   90.0   0.9538570E-05  -0.8978043E-05   0.1410391E-20   0.3185720E-20
C         -0.9164338E-05   0.8979190E-05   0.5112782E-19   0.8300842E-20
C          0.1199542E-19  -0.1247253E-19  -0.1510833E-07   0.1523647E-07
C          0.3839321E-22  -0.1357579E-20  -0.1524935E-07  -0.1335003E-07
C  180.0   0.1850583E-04   0.0000000E+00   0.1005018E-18   0.1096689E-19
C          0.0000000E+00   0.1776269E-04   0.1005018E-18   0.1096689E-19
C          0.2512546E-19   0.6226571E-22  -0.1776269E-04   0.1530543E-21
C         -0.2741722E-20  -0.2450208E-20  -0.4427613E-21  -0.1850583E-04
C
C-----------------------------------------------------------------------
C
      parameter (nmp=np*(np+2),nmp0=(np+1)*(np+4)/2,np2=2*np)
      parameter (NXMAX=3000,nangmax=1801)
C----------------------------------------------------------------------- 
C  NXMAX - the maximum dimension of an auxiliary array in calculating 
C          Mie scattering coefficients, must not be less than 
C                  1.1*np*|m|+10 
C          where |m| is the largest refractive index of all spheres) 
C-----------------------------------------------------------------------
      parameter (ni0=np*(np+1)*(np2+1)/3+np*np)
      parameter (ng0=np*(2*np**3+10*np**2+19*np+5)/6)	
      parameter (nrc=4*np*(np+1)*(np+2)/3+np)
      parameter (nij=nLp*(nLp-1)/2)
      integer u,v,u0,nmax(nLp),uvmax(nLp),ind(nLp)
      double precision k,r0(6,nLp),x(nLp),dang(nangmax),r00(3,nLp),
     +   rsr0(NXMAX),rsi0(NXMAX),rsx0(NXMAX),px0(NXMAX),w1(np),w2(np),
     +   w3(np),w4(np),mue(4,4,nangmax),besj(0:np2+1),besy(0:np2+1),
     +   i11(nangmax),i21(nangmax),i22(nangmax),i12(nangmax),
     +   inat(nangmax),pol(nangmax),cexti(nLp),cabsi(nLp),cscai(nLp),
     +   assymi(nLp),cpri(nLp),drot(nrc,nij),c0i(nLp),c1i(nLp),
     +   bes0(nij),confg(5,nij),taup(2),taupj(2),taupg(2),taupjg(2),
     +   tau0p(2),tau1p(2),tau2p(2),tau0pg(2),tau1pg(2),tau2pg(2),
     +   tau0pj(2),tau1pj(2),tau2pj(2),tau0pjg(2),tau1pjg(2),
     +   tau2pjg(2),tau20(2,2),tau11(2,2),tau02(2,2),tau20g(2,2),
     +   tau11g(2,2),tau02g(2,2),taum(2,2),taumg(2,2),w01s(np2+1),
     +   wmf1(np,np,0:np,0:np2),wm1(np,np,0:np,0:np2),wcf(np2+1),
     +   wsdt(np,np,-np:np,0:np2)
      complex*16 A,B,cmz,Aj,Bj,A2,B2,Aj2,Bj2,A0,B0,ephi,ci,cin,ci0,
     +   A1,B1,Aj1,Bj1,Aj0,Bj0,cmzj,cmzg,Ag,Bg,
     +   A0p(2,2,2,2),A1p(2,2,2,2),B0p(2,2,2,2),B1p(2,2,2,2),
     +   A0pg(2,2,2,2),A1pg(2,2,2,2),B0pg(2,2,2,2),B1pg(2,2,2,2),
     +   atr0(ni0,nij),btr0(ni0,nij),atr(2,np,nmp),at(nmp),bt(nmp),
     +   atr1(ni0,nij),btr1(ni0,nij),ek(np,nij),ref(nLp),ref0(nLp),
     +   p0(nLp,nmp),q0(nLp,nmp),an(np),bn(np),aMie(nLp,np),
     +   bMie(nLp,np),B2i(nLp),at0(nLp,nmp),bt0(nLp,nmp),at1(nmp),
     +   bt1(nmp),as(nLp,nmp),bs(nLp,nmp),as2(nLp,nmp),bs2(nLp,nmp),
     +   tta(nLp,nLp,nmp),ttb(nLp,nLp,nmp),tta0(nLp,nLp,nmp),
     +   ttb0(nLp,nLp,nmp),asr(nLp,nLp,nmp),bsr(nLp,nLp,nmp),
     +   as0(nLp,nLp,nmp),bs0(nLp,nLp,nmp),asc(nLp,nLp,nmp),
     +   bsc(nLp,nLp,nmp),as1(nLp,nLp,nmp),bs1(nLp,nLp,nmp),
     +   ast(nLp,nLp,nmp),bst(nLp,nLp,nmp),asp(nLp,nLp,nmp),
     +   bsp(nLp,nLp,nmp),asv(nLp,nLp,nmp),bsv(nLp,nLp,nmp),
     +	 atj(nmp),btj(nmp),pct(nLp,nLp,2,2,nmp,nmp),
     +   A1m(0:np,np,np,2,2),A2m(0:np,np,np,2,2),
     +   B11n(np,np,2,2),B12n(np,np,2,2),B13n(np,np,2,2),
     +   B21n(np,np,2,2),B22n(np,np,2,2),B23n(np,np,2,2),
     +   fhmf1(0:np2,np,0:np,2),fmf1(0:np2,np,0:np,2),
     +   fhm1v(0:np2,np,0:np,2),fm1v(0:np2,np,0:np,2),
     +   fhm1q(0:np2,np,0:np,2),fm1q(0:np2,np,0:np,2),
     +   fhmf1vq(0:np2,np,0:np,2),fmf1vq(0:np2,np,0:np,2),
     +   fhas(0:np2,2,2),fnhs(0:np2,2,2)
      CHARACTER FLNAME*20,fileout*20,fileout1*19,
     +   tailn*3,cnr2*2,cnr3*3,cnr1*1,flout*22
	COMMON/MIESUB/ twopi,pih
        common/rot/bcof(0:np+2),dc(-np:np,0:nmp)
        common/fnr/fnr(0:4*(np+1))
	common/pitau/pi(nmp0),tau(nmp0)
        common/tran/atr
	common/ig0/iga0(ni0)
	common/g0/ga0(ng0)
	common/cofmnv0/cof0(ni0)
	common/crot/cofsr(nmp)
	pih   = dacos(0.d0)
	twopi  = 4.d0*pih
	pione  = 2.d0*pih
	ci=dcmplx(0.d0,1.d0)
	cin=dcmplx(0.d0,-1.d0)
	ci0=dcmplx(0.d0,0.d0)
	gcs=0.d0
	gcv=0.d0
	eps=1.d-20
	fint=0.01d0
C------------------------------------------------------------------------
C  fint: the interaction index in the range of [0,1] (normally fint=0.01)
C  In the scattering calculations, a quantity "f" for a pair of component 
C  spheres is defined by f=(r_i+r_j)/d_{ij}, where (r_i,r_j) are the radii 
C  of the two spheres and d_{ij} is the separation distance of the two 
C  sphere-centers. When f<fint, interaction between the pair is considered
C  to be negligible and no interaction will be calculated between the two
C  spheres. When fint=0, no sphere is excluded in interaction calculations.
C------------------------------------------------------------------------
	OPEN(UNIT=1,FILE='gmmTr.in',STATUS='OLD')
	READ(1,'(a20)') FLNAME
C-----------------------------------------------------------------------
C  FLNAME - the input file name for the aggregate configuration 
C           the first line in this file is the incident wavelength, 
C           the second line is the number of spheres, 
C           rest lines provide coordinates of sphere-centers, radii, and 
C           refractive indexes of the spheres; each line contains six 
C           real numbers of x, y, z, r, Re(m), Im(m) for a sphere
C-----------------------------------------------------------------------
	write(6,'(a12,a20)') 'input file: ',FLNAME    
        READ(1,*) idMie
C------------------------------------------------------------------------
C  default: idMie=0
C  idMie=1: calculating only coherent Mie-scattering, no interaction 
C------------------------------------------------------------------------
        write(6,'(a7,i3)') 'idMie: ',idMie
   	READ(1,*) small,MXINT
C------------------------------------------------------------------------ 
C  small: error tolerance for the iterative solution process for solving 
C         the interacting scattering coefficients (1.d-6)
C  MXINT: the maximum number of iterations allowed in the iterative 
C  solution process
C  This code uses BI-CGSTAB, the stabilized Bi-Conjugate Gradient method
C  [see H.A. van der Vorst, SIAM J. Sci. Stat. Comput. 13, 631, (1992); 
C  M.H. Gutknecht, SIAM J. Sci. comput. 14, pp.1020-1033 (1993)].
C------------------------------------------------------------------------
	write(6,'(a22,e10.2)') 'Convergence criterion:',small
	write(6,'(a37,i5)') 
     +     'Maximum number of iteration 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]   (Normally, set NADD=0)
C------------------------------------------------------------------------
	write(6,'(a41,i3)') 
     +     'Scat. orders added to Wiscombe criterion:', NADD        
	READ(1,*) idscmt,sang
C------------------------------------------------------------------------ 
C  idscmt: switch for calculation of scattering matrix
C          (idscmt<0 for not calculating scattering matrix)
C  sang: the scattering angle interval for output  
C  example: when sang=1, the results in output will be written for every 
C  degree of scattering angle from 0 to 180 
C------------------------------------------------------------------------
        if(sang.lt.1.d-3) sang=90.d0
	nang=90.d0/sang+1.d0
	nang2=2*nang-1
	if(nang2.gt.nangmax) then
	   write(6,*) 'sang too small'
	   write(6,*) 'please increase sang in the input file gmm01f.in'
	   write(6,*) '  and try again, or'
	   write(6,*) '  increase nangmax in the parameter line of the'
	   write(6,*) '  main code, recompile, then try again'
	   stop
	endif
	close(1)
	write(6,'(/)')	
	OPEN(UNIT=2,FILE=FLNAME,STATUS='OLD')
	READ(2,*) w
C------------------------------------------------------------------------
C  w -- incident wavelength
C------------------------------------------------------------------------
	READ(2,*) nL
C------------------------------------------------------------------------
C  nL -- number of spheres in the aggregate
C------------------------------------------------------------------------
	if(nL.gt.nLp) then
	  write(6,*) 'Parameter nLp too small, must be >', nL 
	  write(6,*) 'Change nLp in gmm01f.par, recompile, then try again'
	  stop
	endif
	if(nL.eq.1) then 
	   idMie=1
	endif
C------------------------------------------------------------------------
C  input the configuration and particle parameters for the aggregate
C  each line includes 6 numbers:
C  x-, y-, z-coordinates of the sphere-center, the radius of the sphere 
C  in the same unit of the incident wavelength, the real and imaginary 
C  parts of the refractive index 
C------------------------------------------------------------------------
	do 1 i=1,nL
	   read(2,*,err=10) (r0(j,i),j=1,6)
	   x0=r0(1,i)
	   y0=r0(2,i)
	   z0=r0(3,i)
	   r00(1,i)=x0
	   r00(2,i)=y0
	   r00(3,i)=z0
	   if(r0(6,i).gt.0.d0) r0(6,i)=-r0(6,i)
	   if(r0(5,i).eq.1.d0.and.r0(6,i).eq.0.d0) goto 1 
	   gcs=gcs+r0(4,i)*r0(4,i)
	   gcv=gcv+r0(4,i)*r0(4,i)*r0(4,i)
 1	continue
        close(2)
	gcsr=dsqrt(gcs)
	gcvr=gcv**(1.d0/3.d0)
	goto 11
 10	write(6,*) 'fatal error in the input file'
	stop
 11	k=twopi/w
        xv=k*gcvr
        xs=k*gcsr
        write(6,'(a,f7.3,a,f7.3)') ' volume-equiv. xv: ',xv,
     +      '  surface-equiv. xs: ',xs
        write(6,'(/)')
        
	fileout='gmm01TrA.out'
	fileout1=fileout
	if(idMie.eq.1) then
	   fileout='M'//fileout1
           write(6,*) '*** Calculating only coherent Mie-scattering ***'
           write(6,*) '*** No interaction included ********************'
	endif

	do i=1,nL
	   x(i)=k*r0(4,i)
	   ref(i)=dcmplx(r0(5,i),r0(6,i))
	   ref0(i)=dcmplx(r0(5,i),-r0(6,i))
	enddo
	do j=1,np
	   do i=1,nL
	      aMie(i,j)=0.d0
	      bMie(i,j)=0.d0
	   enddo
	enddo
	nmax0=1
	do i=1,nL
	   if(i.eq.1) goto  12
	   if(x(i).eq.x(i-1).and.ref(i).eq.ref(i-1)) then
	      nmax(i)=nmax(i-1)
	      uvmax(i)=uvmax(i-1)
	      goto 15
	   endif
 12	   write(6,'(a,i3,a,f10.4)') 'sphere #',i, 
     +         '   individual size parameter: ',x(i)
c
c  calculating Mie-scattering coefficients for each spheres
c  the ratio method of Wang and van der Hulst is used in calculating  
c  Riccati-Bessel functions [see Wang and van der Hulst, Appl. Opt. 
c  30, 106 (1991), Xu, Gustafson, Giovane, Blum, and Tehranian, 
c  Phys. Rev. E 60, 2347 (1999)]
c 
 	   call abMiexud(x(i),ref(i),np,NXMAX,nmax(i),an,bn,NADD,
     +                   rsr0,rsi0,rsx0,px0,w1,w2,w3,w4,eps)
	   if(nmax(i).gt.np) then
	      write(6,*) ' Parameter np too small, must be >',nmax(i)
	      write(6,*) ' Please change np in gmm01f.par, recompile,' 
	      write(6,*) '   then try again'
	      stop
	   endif
	   uvmax(i)=nmax(i)*(nmax(i)+2)
	   write(6,'(a,1x,i4)') 
     +        ' Actual single-sphere expansion truncation:',nmax(i)
	   do j=1,nmax(i)
	      temp1=an(j)
	      temp2=bn(j)
	      if(j.eq.1.or.j.eq.nmax(i)) 
     +           write(6,'(i10,4e15.7)') j,temp1,
     +                 dimag(an(j)),temp2,dimag(bn(j))
	   enddo
 15	   do j=1,nmax(i)
	      aMie(i,j)=an(j)
	      bMie(i,j)=bn(j)
	   enddo
	   if(nmax(i).gt.nmax0) nmax0=nmax(i)
	enddo
	iram=0
	write(6,'(/)')
	write(6,*) 'input sphere-positions: '
	i=1
	write(6,'(i5,3f14.5)') i,r0(1,1),r0(2,1),r0(3,i)
	i=nL
	write(6,'(i5,3f14.5)') i,r0(1,i),r0(2,i),r0(3,i)
C------------------------------------------------------------------------
C  calculating constants and Gaunt coefficients
C------------------------------------------------------------------------
      n0=nmax0+2
      fnr(0)=0.d0
      do n=1,4*(nmax0+1)
         fnr(n)=dsqrt(dble(n))
      enddo
      bcof(0)=1.d0
      do n=0,n0-1
         bcof(n+1)=fnr(n+n+2)*fnr(n+n+1)*bcof(n)/fnr(n+1)/fnr(n+1)
      enddo
C
C  the formulation used here for the calculation of Gaunt coefficients 
C  can be found in Bruning and Lo, IEEE Trans. Anten. Prop. Ap-19, 378 
C  (1971) and Xu, J. Comput. Appl. Math. 85, 53 (1997), J. Comput. Phys. 
C  139, 137 (1998)
C
        call cofsrd(nmax0)	
	call cofd0(nmax0)
        call cofnv0(nmax0)
        call gau0(nmax0)
C------------------------------------------------------------------------
C  calculating rotational and translation coefficients
C------------------------------------------------------------------------
 	do i=1,nL-1
           do j=i+1,nL
              ij=(j-1)*(j-2)/2+j-i
              x0=r0(1,i)-r0(1,j)
              y0=r0(2,i)-r0(2,j)
              z0=r0(3,i)-r0(3,j)
              call carsphd(x0,y0,z0,d,xt,sphi,cphi)
              temp=(r0(4,i)+r0(4,j))/d
              if(temp.lt.1.d0) goto 151
              temp1=r0(4,i)
              if(r0(4,j).lt.r0(4,i)) temp1=r0(4,j)
              temp1=(r0(4,i)+r0(4,j)-d)/temp1
              if(temp1.gt.0.001d0) then
                 write(6,'(/)')
                 write(6,*) 'Fatal error in aggregate configuration:'
                 write(6,*) '        SPHERES OVERLAPPED' 
                 write(6,'(e12.4)') temp1
                 write(6,'(i6,4e12.4)') i,(r0(ilist,i),ilist=1,4)
                 write(6,'(i6,4e12.4)') j,(r0(ilist,j),ilist=1,4)
                 write(6,*) 'This code calculates aggregate-scattering'
                 write(6,*) 'of non-intersecting spheres. The ratio of'
                 write(6,*) '(r1+r2)/d must not exceed 1, where r1+r2'
                 write(6,*) 'is the sum of the radii of any pair of'
                 write(6,*) 'spheres and d is the separation distance' 
                 write(6,*) 'between the two sphere centers. This ratio'
                 write(6,*) 'for the pair of spheres listed is greater'
                 write(6,*) 'than 1 and the overlapped portion exceeds'
                 write(6,*) '0.001 of the (smaller) sphere-radius as '
                 write(6,*) 'shown above. The sequence numbers and the'
                 write(6,*) 'Cartesian coordinates of the two sphere' 
                 write(6,*) 'centers, the radii of the two spheres are'
                 write(6,*) 'also shown. Please correct the input data' 
                 write(6,*) 'for the aggregate configuration and then'
                 write(6,*) 'run again.'
                 stop
              endif
 151          confg(1,ij)=x0
              confg(2,ij)=y0
              confg(3,ij)=z0
              confg(4,ij)=d
              confg(5,ij)=temp
              if(temp.le.fint) goto 16            
              ephi=dcmplx(cphi,sphi)
              nlarge=max(nmax(i),nmax(j))
              do m=1,nlarge
                 ek(m,ij)=ephi**m
              enddo
              xd=k*d
              nbes=2*nlarge+1
	      call besseljd(nbes,xd,besj)
	      call besselyd(nbes,xd,besy)
C
C  calculating (-1)^{m+u)d_{mu}^n where d_{mu}^n is the "reduced rotation 
C  matrix elements"  
C
              call rotcoef(xt,nlarge)
              irc=0
              do n=1,nlarge
                 n1=n*(n+1)
                 do u=-n,n
                    do m=-n,n
                       imn=n1+m
                       irc=irc+1
                       drot(irc,ij)=dc(u,imn)
                    enddo
                 enddo
              enddo
              itrc=0
              nsmall=min(nmax(i),nmax(j))
C
C  the formulation used here for the calculation of vector translation 
C  coefficients are from Cruzan, Q. Appl. Math. 20, 33 (1962) and  
C  Xu, J. Comput. Phys. 139, 137 (1998)
C
              do m=-nsmall,nsmall
                 n1=max(1,iabs(m))
                 do n=n1,nlarge
                    do v=n1,nlarge
                       itrc=itrc+1
                       call cofxuds0(nmax0,m,n,v,besj,besy,
     +                         atr0(itrc,ij),btr0(itrc,ij),
     +                         atr1(itrc,ij),btr1(itrc,ij))     
                    enddo
                 enddo
              enddo
 16        continue
           enddo
        enddo
	if(idMie.eq.1) then
	   do j=1,nL
 	      do iuv=1,uvmax(j)
 	         v=dsqrt(dble(iuv))
 	         do i=1,nL
 	            do imn=1,uvmax(i)
	               pct(j,i,1,1,imn,iuv)=0.d0
	               pct(j,i,1,2,imn,iuv)=0.d0
	               pct(j,i,2,1,imn,iuv)=0.d0
	               pct(j,i,2,2,imn,iuv)=0.d0
	               if(j.eq.i.and.imn.eq.iuv) then 
	                  pct(j,i,1,1,imn,iuv)=aMie(j,v)
	                  pct(j,i,2,2,imn,iuv)=bMie(j,v)
	               endif     
 	            enddo
 	         enddo
 	      enddo
 	   enddo
           goto 1800
        endif
C--------------------------------------------------------------------
c  begins iteration process to solve T-matrix columns of (u,v,q)
c  BI-CGSTAB [see van der Vorst, SIAM J. Sci. Stat. Comput. 13, 631 
c  (1992); Gutknecht, SIAM J. Sci. comput. 14, pp. 1020-1033 (1993)]   
C--------------------------------------------------------------------
      do i=1,nL-1
         do j=i+1,nL
            ij=(j-1)*(j-2)/2+j-i
            d=confg(4,ij)
            xd=k*d
            bes0(ij)=dsin(xd)/xd	      
         enddo
      enddo
      write(6,*) 'Starting Bi-CGSTAB to solve T-matrix'         
      n0=nmax0*(nmax0+2)
      do 1001 iuv=1,n0
         v=dsqrt(dble(iuv))
         iuvc=v*v
c	 iuvc=1
         u=iuv-v*(v+1)
         do 1002 iq=1,2
            iuv1=-u+v+v*v
            do imn=1,iuvc-1
               n=dsqrt(dble(imn))
               m=imn-n*n-n
               imn1=-m+n*n+n
               cz=(-1)**(m+u+n+v)
               do j=1,nL
                  do i=1,nL
                     asr(j,i,imn)=cz*pct(i,j,iq,1,iuv1,imn1)
                     bsr(j,i,imn)=cz*pct(i,j,iq,2,iuv1,imn1)
                  enddo
               enddo
            enddo
            do j=1,nL
               ind(j)=0
               do i=1,nL
                  do imn=1,iuvc-1
                     as(i,imn)=asr(j,i,imn)
                     bs(i,imn)=bsr(j,i,imn)
                  enddo
                  do imn=iuvc,uvmax(i)
                     as(i,imn)=ci0
                     bs(i,imn)=ci0
                  enddo
               enddo
               call transT(nL,r0,nmax,uvmax,fint,atr0,btr0,
     +                 ek,drot,as,bs,as2,bs2,ind,confg,iuvc,2)
               do i=1,nL
                  do imn=iuvc,uvmax(i)
                     tta0(j,i,imn)=as2(i,imn)
                     ttb0(j,i,imn)=bs2(i,imn)
                  enddo
               enddo
            enddo	      	      	      	      
            do j=1,nL
               do i=1,nL
                  do imn=iuvc,n0
                     asr(j,i,imn)=ci0
                     bsr(j,i,imn)=ci0
                  enddo
               enddo
            enddo
            do j=1,nL
               if(iuv.le.uvmax(j)) then
                  if(iq.eq.1) then
                     asr(j,j,iuv)=aMie(j,v)
                  else
                     bsr(j,j,iuv)=bMie(j,v)
                  endif
               endif
            enddo
            do j=1,nL
               if(iq.eq.1) then
                  c0i(j)=aMie(j,v)*dconjg(aMie(j,v))
               else
                  c0i(j)=bMie(j,v)*dconjg(bMie(j,v))
               endif
            enddo
            niter=1
            do j=1,nL
               do i=1,nL
                  do imn=1,iuvc-1
                     as(i,imn)=ci0
                     bs(i,imn)=ci0
                  enddo
                  do imn=iuvc,uvmax(i)
                     as(i,imn)=asr(j,i,imn)
                     bs(i,imn)=bsr(j,i,imn)
                  enddo
               enddo
               call transT(nL,r0,nmax,uvmax,fint,atr0,btr0,
     +              ek,drot,as,bs,as2,bs2,ind,confg,iuvc,1)
               do i=1,nL
                  do imn=iuvc,uvmax(i)
                     tta(j,i,imn)=as2(i,imn)+tta0(j,i,imn)
                     ttb(j,i,imn)=bs2(i,imn)+ttb0(j,i,imn)
                  enddo
               enddo
            enddo
            do 611 i=1,nL
               c1i(i)=0.d0
               do 6111 imn=iuvc,uvmax(i)
                  n=dsqrt(dble(imn))
                  do 6112 j=1,nL
                     if(iuv.gt.uvmax(j)) goto 6112
                     as1(j,i,imn)=0.d0
                     bs1(j,i,imn)=0.d0
                     if(imn.eq.iuv) then
                        ijmax=max(i,j)
                        ijmin=min(i,j)
                        if(ijmax.eq.ijmin) then
                           cz=1.d0
                        else
                           ij=(ijmax-1)*(ijmax-2)/2+ijmax-ijmin
                           cz=bes0(ij)
                        endif
                        if(iq.eq.1) then
                           as1(j,i,imn)=as1(j,i,imn)+cz*aMie(i,n)
                        else
                           bs1(j,i,imn)=bs1(j,i,imn)+cz*bMie(i,n)
                        endif
                     endif
                     do 6113 jj=1,nL
                        ijmax=max(jj,j)
                        ijmin=min(jj,j)
                        if(ijmax.eq.ijmin) then
                           cz=1.d0
                        else
                           ij=(ijmax-1)*(ijmax-2)/2+ijmax-ijmin
                           cz=bes0(ij)
                        endif
                        A0=asr(jj,i,imn)
                        B0=bsr(jj,i,imn)
                        A0=A0+aMie(i,n)*tta(jj,i,imn)
                        B0=B0+bMie(i,n)*ttb(jj,i,imn)
                        as1(j,i,imn)=as1(j,i,imn)-cz*A0
                        bs1(j,i,imn)=bs1(j,i,imn)-cz*B0
 6113                continue
                     A=as1(j,i,imn)
                     B=bs1(j,i,imn)
                     c1i(i)=c1i(i)+A*dconjg(A)
                     c1i(i)=c1i(i)+B*dconjg(B)
 6112	          continue
 6111 	       continue
 611        continue
            temp=0.d0
            B0=0.d0
            do 612 i=1,nL
               cext0=c1i(i)/c0i(i)
               if(cext0.lt.small) ind(i)=1
               if(ind(i).gt.0) goto 612
               if(cext0.gt.temp) temp=cext0
               B0=B0+c1i(i)
 612	    continue
            if(temp.lt.small) then
               do i=1,nL
                  do imn=1,uvmax(i)
                     do j=1,nL
                        pct(j,i,1,iq,imn,iuv)=asr(j,i,imn)
                        pct(j,i,2,iq,imn,iuv)=bsr(j,i,imn)
                     enddo
                  enddo
               enddo
               goto 1002
            endif
            do 613 i=1,nL
               if(ind(i).gt.0) goto 613
               do imn=1,iuvc-1
                  do j=1,nL
                     asp(j,i,imn)=ci0
                     bsp(j,i,imn)=ci0
                     as0(j,i,imn)=ci0
                     bs0(j,i,imn)=ci0
                  enddo
               enddo
               do imn=iuvc,uvmax(i)
                  do j=1,nL
                     asp(j,i,imn)=as1(j,i,imn)
                     bsp(j,i,imn)=bs1(j,i,imn)
                     as0(j,i,imn)=as1(j,i,imn)
                     bs0(j,i,imn)=bs1(j,i,imn)
                  enddo
               enddo
 613        continue
            do j=1,nL
               do i=1,nL
                  do imn=1,iuvc-1
                     as(i,imn)=ci0
                     bs(i,imn)=ci0
                  enddo
                  do imn=iuvc,uvmax(i)
                     as(i,imn)=asp(j,i,imn)
                     bs(i,imn)=bsp(j,i,imn)
                  enddo
               enddo
               call transT(nL,r0,nmax,uvmax,fint,atr0,btr0,
     +              ek,drot,as,bs,as2,bs2,ind,confg,iuvc,1)
               do i=1,nL
                  do imn=iuvc,uvmax(i)
                     tta(j,i,imn)=as2(i,imn)
                     ttb(j,i,imn)=bs2(i,imn)
                  enddo
               enddo
            enddo
            A0=0.d0
            do 614 i=1,nL
               if(ind(i).gt.0) goto 614
               do 6141 imn=iuvc,uvmax(i)
                  n=dsqrt(dble(imn))
                  do 6142 j=1,nL
                     if(iuv.gt.uvmax(j)) goto 6142
                     ast(j,i,imn)=0.d0
                     bst(j,i,imn)=0.d0
                     do 6143 jj=1,nL
                        ijmax=max(jj,j)
                        ijmin=min(jj,j)
                        if(ijmax.eq.ijmin) then
                           cz=1.d0
                        else
                           ij=(ijmax-1)*(ijmax-2)/2+ijmax-ijmin
                           cz=bes0(ij)
                        endif
                        Aj2=asp(jj,i,imn)
                        Bj2=bsp(jj,i,imn)
                        Aj2=Aj2+aMie(i,n)*tta(jj,i,imn)
                        Bj2=Bj2+bMie(i,n)*ttb(jj,i,imn)
                        ast(j,i,imn)=ast(j,i,imn)+cz*Aj2
                        bst(j,i,imn)=bst(j,i,imn)+cz*Bj2
 6143                continue	          
                     A0=A0+dconjg(as0(j,i,imn))*ast(j,i,imn)
                     A0=A0+dconjg(bs0(j,i,imn))*bst(j,i,imn)
 6142             continue
 6141          continue
 614        continue
            if(cdabs(A0).lt.1.d-200) then
               do i=1,nL
                  do imn=1,uvmax(i)
                     do j=1,nL
                        pct(j,i,1,iq,imn,iuv)=asr(j,i,imn)
                        pct(j,i,2,iq,imn,iuv)=bsr(j,i,imn)
                     enddo
                  enddo
               enddo  
               goto 1002
            endif
            Aj=B0/A0
  62        do 621 i=1,nL
               if(ind(i).gt.0) goto 621
               do imn=iuvc,uvmax(i)
                  do 6211 j=1,nL
                     if(iuv.gt.uvmax(j)) then
                        asv(j,i,imn)=0.d0
                        bsv(j,i,imn)=0.d0
                        go to 6211
                     endif
                     asv(j,i,imn)=asp(j,i,imn)-Aj*ast(j,i,imn)
                     bsv(j,i,imn)=bsp(j,i,imn)-Aj*bst(j,i,imn)
 6211             continue 
               enddo
 621        continue
            do j=1,nL
               do i=1,nL
                  do imn=1,iuvc-1
                     as(i,imn)=ci0
                     bs(i,imn)=ci0
                  enddo
                  do imn=iuvc,uvmax(i)
                     as(i,imn)=asv(j,i,imn)
                     bs(i,imn)=bsv(j,i,imn)
                  enddo
               enddo
               call transT(nL,r0,nmax,uvmax,fint,atr0,btr0,
     +              ek,drot,as,bs,as2,bs2,ind,confg,iuvc,1)
               do i=1,nL
                  do imn=iuvc,uvmax(i)
                     tta(j,i,imn)=as2(i,imn)
                     ttb(j,i,imn)=bs2(i,imn)
                  enddo
               enddo
            enddo
            A2=0.d0
            B2=0.d0
            do 622 i=1,nL
               if(ind(i).gt.0) goto 622
               do 6221 imn=iuvc,uvmax(i)
                  n=dsqrt(dble(imn))
                  do 6222 j=1,nL
                     if(iuv.gt.uvmax(j)) goto 6222
                     asc(j,i,imn)=0.d0
                     bsc(j,i,imn)=0.d0
                     do 6223 jj=1,nL
                        ijmax=max(jj,j)
                        ijmin=min(jj,j)
                        if(ijmax.eq.ijmin) then
                           cz=1.d0
                        else
                           ij=(ijmax-1)*(ijmax-2)/2+ijmax-ijmin
                           cz=bes0(ij)
                        endif
                        Aj2=asv(jj,i,imn)
                        Bj2=bsv(jj,i,imn)
                        Aj2=Aj2+aMie(i,n)*tta(jj,i,imn)
                        Bj2=Bj2+bMie(i,n)*ttb(jj,i,imn)
                        asc(j,i,imn)=asc(j,i,imn)+cz*Aj2
                        bsc(j,i,imn)=bsc(j,i,imn)+cz*Bj2
 6223                continue	            
                     A2=A2+dconjg(asc(j,i,imn))*asv(j,i,imn)
                     A2=A2+dconjg(bsc(j,i,imn))*bsv(j,i,imn)
                     B2=B2+dconjg(asc(j,i,imn))*asc(j,i,imn)
                     B2=B2+dconjg(bsc(j,i,imn))*bsc(j,i,imn)
 6222             continue
 6221          continue
 622        continue
            if(cdabs(B2).lt.1.d-200) then
               do i=1,nL
                  do imn=1,uvmax(i)
                     do j=1,nL
                        pct(j,i,1,iq,imn,iuv)=asr(j,i,imn)
                        pct(j,i,2,iq,imn,iuv)=bsr(j,i,imn)
                     enddo
                  enddo
               enddo  
               goto 1002
            endif
            Bj=A2/B2
            do 623 i=1,nL
               if(ind(i).gt.0) goto 623
               do imn=iuvc,uvmax(i)
                  do 6231 j=1,nL
                     if(iuv.gt.uvmax(j)) then
                        asp(j,i,imn)=ci0
                        bsp(j,i,imn)=ci0
                        go to 6231
                     endif   
                     asp(j,i,imn)=asv(j,i,imn)-Bj*asc(j,i,imn)
                     bsp(j,i,imn)=bsv(j,i,imn)-Bj*bsc(j,i,imn)
 6231             continue
               enddo
 623        continue
            do 624 i=1,nL
               c1i(i)=0.d0
               if(ind(i).gt.0) goto 624
               do 6241 imn=iuvc,uvmax(i)
                  do 6242 j=1,nL
                     if(iuv.gt.uvmax(j)) goto 6242
                     Aj2=Aj*as1(j,i,imn)+Bj*asv(j,i,imn)
                     Bj2=Aj*bs1(j,i,imn)+Bj*bsv(j,i,imn)
                     asr(j,i,imn)=asr(j,i,imn)+Aj2
                     bsr(j,i,imn)=bsr(j,i,imn)+Bj2
                     c1i(i)=c1i(i)+Aj2*dconjg(Aj2)
                     c1i(i)=c1i(i)+Bj2*dconjg(Bj2)
 6242             continue
 6241          continue
 624        continue
            cext0=0.d0
            cext1=0.d0
            do 625 i=1,nL
               if(ind(i).gt.0) goto 625
               cext0=cext0+c0i(i)
               cext1=cext1+c1i(i)
 625        continue
            temp=cext1/cext0
            if(iuv.eq.1.and.iq.eq.1) then
               if(niter/20*20.eq.niter.or.niter.eq.1) then
                  write(6,'(a11,i4,2x,e15.7)')
     +                 'iteration #',niter,temp
               endif
            endif
            if(temp.lt.small) then
               if(u.eq.v) then
                  write(6,'(3i4,a12,i4,2x,e15.5)') u,v,iq,
     +                 '  iteration ',niter,temp	         
               endif
               do i=1,nL
                  do imn=1,uvmax(i)
                     do j=1,nL
                        pct(j,i,1,iq,imn,iuv)=asr(j,i,imn)
                        pct(j,i,2,iq,imn,iuv)=bsr(j,i,imn)
                     enddo
                  enddo
               enddo  
               goto 1002
            endif
            if(niter.gt.MXINT) then
               write(6,'(a9,i5,i4,i3)') 'Caution: ',u,v,iq
               write(6,*) '*** Maximum iterations exceeded ***'
               write(6,*) '*** Solution may be inaccurate  ***'	
               goto 1002
            endif
            B2=0.d0
            do 626 i=1,nL 	   
               if(ind(i).gt.0) goto 626
               B2i(i)=0.d0
               do imn=iuvc,uvmax(i)
                  do 6261 j=1,nL
                     if(iuv.gt.uvmax(j)) goto 6261
                     Aj2=dconjg(as0(j,i,imn))*asp(j,i,imn)
                     Bj2=dconjg(bs0(j,i,imn))*bsp(j,i,imn)
                     B2i(i)=B2i(i)+Aj2
                     B2i(i)=B2i(i)+Bj2
 6261             continue
               enddo
 	       B2=B2+B2i(i)
 626        continue	
            A0=B0*Bj
            if(cdabs(A0).lt.1.d-200) then
               do i=1,nL
                  do imn=1,uvmax(i)
                     do j=1,nL
                        pct(j,i,1,iq,imn,iuv)=asr(j,i,imn)
                        pct(j,i,2,iq,imn,iuv)=bsr(j,i,imn)
                     enddo
                  enddo
               enddo  
               goto 1002
            endif
            A0=-Aj*B2/A0
            do 627 i=1,nL
               if(ind(i).gt.0) goto 627	   
               do imn=iuvc,uvmax(i)
                  do 6271 j=1,nL
                     if(iuv.gt.uvmax(j)) goto 6271
                     as1(j,i,imn)=as1(j,i,imn)-Bj*ast(j,i,imn)
                     bs1(j,i,imn)=bs1(j,i,imn)-Bj*bst(j,i,imn)
                     as1(j,i,imn)=asp(j,i,imn)-A0*as1(j,i,imn)
                     bs1(j,i,imn)=bsp(j,i,imn)-A0*bs1(j,i,imn)
 6271             continue
               enddo
 627        continue
            B0=B2
            do j=1,nL
               do i=1,nL
                  do imn=1,iuvc-1
                     as(i,imn)=ci0
                     bs(i,imn)=ci0
                  enddo
                  do imn=iuvc,uvmax(i)
                     as(i,imn)=as1(j,i,imn)
                     bs(i,imn)=bs1(j,i,imn)
                  enddo
               enddo
               call transT(nL,r0,nmax,uvmax,fint,atr0,btr0,
     +              ek,drot,as,bs,as2,bs2,ind,confg,iuvc,1)
               do i=1,nL
                  do imn=iuvc,uvmax(i)
                     tta(j,i,imn)=as2(i,imn)
                     ttb(j,i,imn)=bs2(i,imn)
                  enddo
               enddo
            enddo
            A0=0.d0
            do 629 i=1,nL
               if(ind(i).gt.0) goto 629
               do 6291 imn=iuvc,uvmax(i)
                  n=dsqrt(dble(imn))
                  do 6292 j=1,nL
                     if(iuv.gt.uvmax(j)) goto 6292
                     ast(j,i,imn)=0.d0
                     bst(j,i,imn)=0.d0
                     do 6293 jj=1,nL
                        ijmax=max(jj,j)
                        ijmin=min(jj,j)
                        if(ijmax.eq.ijmin) then
                           cz=1.d0
                        else
                           ij=(ijmax-1)*(ijmax-2)/2+ijmax-ijmin
                           cz=bes0(ij)
                        endif
                        Aj2=as1(jj,i,imn)
                        Bj2=bs1(jj,i,imn)
                        Aj2=Aj2+aMie(i,n)*tta(jj,i,imn)
                        Bj2=Bj2+bMie(i,n)*ttb(jj,i,imn)
                        ast(j,i,imn)=ast(j,i,imn)+cz*Aj2
                        bst(j,i,imn)=bst(j,i,imn)+cz*Bj2
 6293                continue	 
                     A0=A0+dconjg(as0(j,i,imn))*ast(j,i,imn)
                     A0=A0+dconjg(bs0(j,i,imn))*bst(j,i,imn)
 6292             continue
 6291          continue
 629        continue
            if(cdabs(A0).lt.1.d-200) then
               do i=1,nL
                  do imn=1,uvmax(i)
                     do j=1,nL
                        pct(j,i,1,iq,imn,iuv)=asr(j,i,imn)
                        pct(j,i,2,iq,imn,iuv)=bsr(j,i,imn)
                     enddo
                  enddo
               enddo  
               goto 1002
            endif
            Aj=B0/A0
            niter=niter+1
            goto 62
 1002    continue
 1001 continue	
C------------------------------------------------------------------------
C  calculating random-orientation averaged total and individual-particle
C  extinction cross-sections
C------------------------------------------------------------------------
 1800 do i=1,nL
         ind(i)=0
         cexti(i)=0.d0
         cscai(i)=0.d0
         cpri(i)=0.d0
      enddo
      cext=0.d0
      csca=0.d0
      cpr=0.d0
      n0=nmax0*(nmax0+2)
      do iuv=1,n0
         do iq=1,2
            do 1801 j=1,nL
               if(iuv.gt.uvmax(j)) goto 1801
               do i=1,nL
                  do imn=1,uvmax(i)
                     as(i,imn)=pct(j,i,1,iq,imn,iuv)
                     bs(i,imn)=pct(j,i,2,iq,imn,iuv)
                  enddo
               enddo
               call transT(nL,r0,nmax,uvmax,fint,atr1,btr1,
     +                     ek,drot,as,bs,as2,bs2,ind,confg,1,1)
               do i=1,nL
                  do imn=1,uvmax(i)
                     A=as2(i,imn)+as(i,imn)
                     B=bs2(i,imn)+bs(i,imn)
                     pct(j,i,1,iq,imn,iuv)=A
                     pct(j,i,2,iq,imn,iuv)=B
                  enddo
               enddo
               cz=pct(j,j,iq,iq,iuv,iuv)
               cext=cext+cz
               cexti(j)=cexti(j)+cz
 1801	    continue
         enddo
      enddo
C------------------------------------------------------------------------
C  calculating random-orientation averaged asymmetry parameter and  
C  total and individual-particle scattering cross-sections 
C------------------------------------------------------------------------
        do j=1,nL
 	   do iuv=1,uvmax(j)
 	      v=dsqrt(dble(iuv))
 	      u=iuv-v*v-v
 	      iuv1=iuv-2*u
 	      cb=(-1)**(v+u)
 	      jv1=v+1
 	      iuv2=-u+jv1*(jv1+1)
 	      iuv3=-u+(v-1)*v
 	      fuv1=v*jv1
 	      fuv1=dble(u)/fuv1
 	      fuv2=fnr(v)*fnr(v+2)/fnr(2*v+1)
	      fuv2=fuv2/fnr(2*v+3)/dble(jv1)
	      fuv2=fuv2*fnr(jv1-u)*fnr(jv1+u)
	      fuv3=fnr(v-1)*fnr(jv1)/fnr(2*v-1)
	      fuv3=fuv3/fnr(2*v+1)/dble(v)
	      fuv3=fuv3*fnr(v-u)*fnr(v+u)
	      juv1=iuv1-1
	      juv2=iuv1+1
	      juv3=iuv2-1
	      juv4=iuv3-1
	      juv5=iuv2+1
	      juv6=iuv3+1
	      guv=v*(v+1)
	      guv1=-fnr(v-u)*fnr(v+u+1)/guv
	      guv2=-fnr(v+u)*fnr(v-u+1)/guv
	      guv=dble(v+1)*fnr(v)*fnr(v+2)*fnr(2*v+1)*fnr(2*v+3)
	      guv3=u*(u+1)+((v-u)*(v+u+3)+(v+u)*(v-u+1))/2
	      guv3=guv3*fnr(v+u+1)*fnr(v+u+2)/guv
	      guv5=u*(u-1)+((v-u)*(v+u+1)+(v+u)*(v-u+3))/2
	      guv5=-guv5*fnr(v-u+1)*fnr(v-u+2)/guv
	      if(v.gt.1) then 
                 guv=dble(v)*fnr(v-1)*fnr(v+1)*fnr(2*v-1)*fnr(2*v+1)
                 guv4=u*(u+1)+((v-u-2)*(v+u+1)+(v+u)*(v-u+1))/2
                 ntemp=v-u-1
                 if(ntemp.lt.0) then
                    temp=0.d0
                 else
                    temp=fnr(ntemp)
                 endif
	         guv4=-guv4*fnr(v-u)*temp/guv
	         guv6=u*(u-1)+((v-u)*(v+u+1)+(v+u-2)*(v-u+1))/2
	         ntemp=v+u-1
                 if(ntemp.lt.0) then
                    temp=0.d0
                 else
                    temp=fnr(ntemp)
                 endif
	         guv6=guv6*fnr(v+u)*temp/guv
	      endif
 	      do i=1,nL
 	         do imn=1,uvmax(i)
 	            n=dsqrt(dble(imn))
	            m=imn-n*n-n
	            imn1=-m+n*n+n
	            sb=(-1)**(m+n)
	            sb=cb*sb
	            n1=n+1
	            n2=2*n
	            rn=1.0d0/dble(n*n1)
	            p=fnr(n)*fnr(n+2)/fnr(n2+1)
	            p=p/fnr(n2+3)/dble(n1)
	            t=fnr(n-1)*fnr(n+1)/fnr(n2-1)
	            t=t/fnr(n2+1)/dble(n)
	            rm=dble(m)*rn
	            imn2=(n+1)*(n+2)+m
	            fnp=fnr(n+m+1)*fnr(n-m+1)*p	
	            imn3=(n-1)*n+m
 	            fn=fnr(n+m)*fnr(n-m)*t
 	            jmn1=imn+1
 	            jmn2=imn-1
 	            jmn3=imn2+1
 	            jmn4=imn3+1
 	            jmn5=imn2-1
 	            jmn6=imn3-1
 	            gmn=n*(n+1)
	            gmn1=-fnr(n-m)*fnr(n+m+1)/gmn
	            gmn2=-fnr(n+m)*fnr(n-m+1)/gmn
	            gmn=dble(n+1)*fnr(n)*fnr(n+2)
	            gmn=gmn*fnr(2*n+1)*fnr(2*n+3)
	            gmn3=m*(m+1)+((n-m)*(n+m+3)+(n+m)*(n-m+1))/2
	            gmn3=gmn3*fnr(n+m+1)*fnr(n+m+2)/gmn
	            gmn5=m*(m-1)+((n-m)*(n+m+1)+(n+m)*(n-m+3))/2
	            gmn5=-gmn5*fnr(n-m+1)*fnr(n-m+2)/gmn
	            if(n.gt.1) then
	               gmn=dble(n)*fnr(n-1)*fnr(n+1)
	               gmn=gmn*fnr(2*n-1)*fnr(2*n+1)
	               gmn4=m*(m+1)+((n-m-2)*(n+m+1)+(n+m)*(n-m+1))/2
	               ntemp=n-m-1
	               if(ntemp.lt.0) then
                          temp=0.d0
                       else
                          temp=fnr(ntemp)
                       endif
	               gmn4=-gmn4*fnr(n-m)*temp/gmn
	               gmn6=m*(m-1)+((n-m)*(n+m+1)+(n+m-2)*(n-m+1))/2
	               ntemp=n+m-1
	               if(ntemp.lt.0) then
                          temp=0.d0
                       else
                          temp=fnr(ntemp)
                       endif
	               gmn6=gmn6*fnr(n+m)*temp/gmn
	            endif
	            do ip=1,2
	               do iq=1,2
	                  A=pct(j,i,ip,iq,imn,iuv)
	                  B=sb*pct(i,j,iq,ip,iuv1,imn1)
	                  B=dconjg(B)
	                  cz=B*A
	                  csca=csca+cz
	                  cscai(j)=cscai(j)+cz	                  
	                  A=rm*pct(j,i,3-ip,iq,imn,iuv)
	                  if(n.eq.nmax(i)) goto 1951
 	                  A=A+fnp*pct(j,i,ip,iq,imn2,iuv)
 1951	                  if(n.eq.1.or.iabs(m).gt.n-1) goto 1952
 	                  A=A+fn*pct(j,i,ip,iq,imn3,iuv)
 1952                     B=fuv1*pct(i,j,3-iq,ip,iuv1,imn1)
                          if(v.eq.nmax(j)) goto 1953
                          B=B-fuv2*pct(i,j,iq,ip,iuv2,imn1)
 1953                     if(v.eq.1.or.iabs(u).gt.v-1) goto 1954
                          B=B-fuv3*pct(i,j,iq,ip,iuv3,imn1)
 1954	                  B=sb*dconjg(B)
                          cz=B*A
                          cpr=cpr+cz
                          cpri(j)=cpri(j)+cz
                          A=ci0
                          if(iabs(m+1).gt.n) goto 1961
                          A=gmn1*pct(j,i,3-ip,iq,jmn1,iuv)
 1961                     if(n.eq.nmax(i)) goto 1962
                          A=A+gmn3*pct(j,i,ip,iq,jmn3,iuv)
 1962                     if(n.eq.1) goto 1963
                          if(iabs(m+1).gt.n-1) goto 1963
                          A=A+gmn4*pct(j,i,ip,iq,jmn4,iuv)
 1963                     B=ci0
                          if(iabs(u-1).gt.v) goto 1964
                          B=-guv2*pct(i,j,3-iq,ip,juv2,imn1)
 1964                     if(v.eq.nmax(j)) goto 1965
                          B=B+guv5*pct(i,j,iq,ip,juv5,imn1)
 1965                     if(v.eq.1) goto 1966
                          if(iabs(u-1).gt.v-1) goto 1966
                          B=B+guv6*pct(i,j,iq,ip,juv6,imn1)
 1966                     cz=dconjg(B)*A
                          temp=0.5d0*sb*cz
                          cpr=cpr+temp
                          cpri(j)=cpri(j)+temp
                          A=ci0
                          if(iabs(m-1).gt.n) goto 1971
                          A=gmn2*pct(j,i,3-ip,iq,jmn2,iuv)
 1971                     if(n.eq.nmax(i)) goto 1972
                          A=A+gmn5*pct(j,i,ip,iq,jmn5,iuv)
 1972                     if(n.eq.1) goto 1973
                          if(iabs(m-1).gt.n-1) goto 1973
                          A=A+gmn6*pct(j,i,ip,iq,jmn6,iuv)
 1973                     B=ci0
                          if(iabs(u+1).gt.v) goto 1974
                          B=-guv1*pct(i,j,3-iq,ip,juv1,imn1)
 1974                     if(v.eq.nmax(j)) goto 1975
                          B=B+guv3*pct(i,j,iq,ip,juv3,imn1)
 1975                     if(v.eq.1) goto 1976
                          if(iabs(u+1).gt.v-1) goto 1976
                          B=B+guv4*pct(i,j,iq,ip,juv4,imn1)
 1976                     cz=dconjg(B)*A
                          temp=0.5d0*sb*cz
                          cpr=cpr+temp
                          cpri(j)=cpri(j)+temp
	               enddo
	            enddo     
 	         enddo
 	      enddo
 	   enddo
 	enddo 	   
      if(idscmt.lt.0) goto 2000
C------------------------------------------------------------------------
C  calculating random-orientation averaged Mueller matrix elements
C------------------------------------------------------------------------
      write(6,'(/)')
      write(6,*) 'Calculating (Mueller) scattering matrix'
      write(6,'(/)') 
      n0=nmax0*(nmax0+2)
      nmax2=2*nmax0      
      do n=1,nmax0
         do v=1,nmax0
            do m=0,n
               it=n+v+1
               call xuwigd(n,v,m,-1,w01s,wcf,it,nmf1)
               do ii=1,nmf1
                  it=it-1
                  wmf1(n,v,m,it)=w01s(ii)
               enddo
               it=n+v+1
               call xuwigd(n,v,m,1,w01s,wcf,it,nm1)
               do ii=1,nm1
                  it=it-1
                  wm1(n,v,m,it)=w01s(ii)
               enddo
            enddo
         enddo
      enddo
      do n=1,nmax0
         do jn=1,nmax0
            do ip=1,2
               do jp=1,2
                  B11n(n,jn,ip,jp)=ci0
                  B12n(n,jn,ip,jp)=ci0
                  B13n(n,jn,ip,jp)=ci0
                  B21n(n,jn,ip,jp)=ci0
                  B22n(n,jn,ip,jp)=ci0
                  B23n(n,jn,ip,jp)=ci0
                  do m=0,nmax0
                     A1m(m,n,jn,ip,jp)=ci0
                     A2m(m,n,jn,ip,jp)=ci0
                  enddo
               enddo
            enddo
         enddo
      enddo
      do 1900 ids=-nmax2,nmax2         
         do n=1,nmax0          
            do 19001 v=1,nmax0
               if(iabs(ids).gt.n+v) goto 19001
               fv0=0.5d0*fnr(2*v+1)
               do 19002 ms=-n,n
                  mw=ms+ids
                  if(iabs(mw).gt.v) goto 19002
                  it=n+v+1
                  call xuwigd(n,v,ms,-mw,w01s,wcf,it,nvs)
                  fv=((-1)**ms)*fv0
                  do ii=1,nvs
                     it=it-1
                     wsdt(n,v,ms,it)=fv*w01s(ii)
                  enddo
19002          continue
19001       continue
         enddo
         do i=1,nL
            do j=1,nL
               do it=0,nmax2
                  do n=1,nmax0
                     do m=0,n
                        do ip=1,2
                           fhmf1(it,n,m,ip)=ci0
                           fmf1(it,n,m,ip)=ci0
                           fhm1v(it,n,m,ip)=ci0
                           fm1v(it,n,m,ip)=ci0
                           fhm1q(it,n,m,ip)=ci0
                           fm1q(it,n,m,ip)=ci0
                           fhmf1vq(it,n,m,ip)=ci0
                           fmf1vq(it,n,m,ip)=ci0
                        enddo
                     enddo
                  enddo
               enddo
               do 19003 n=1,nmax0
                  if(n.gt.nmax(i)) goto 19003
                  in0=n*(n+1)
                  cv=1.d0
                  do 19004 v=1,nmax0
                     if(v.gt.nmax(j)) goto 19004
                     if(iabs(ids).gt.n+v) goto 19004
                     cv=-cv
                     iv0=v*(v+1)
                     do it=0,n+v
                        do ip=1,2
                           do iq=1,2
                              fhas(it,ip,iq)=ci0
                              fnhs(it,ip,iq)=ci0
                           enddo
                        enddo
                     enddo
                     do 1901 ms=-n,n
                        mw=ms+ids
                        if(iabs(mw).gt.v) goto 1901
                        it=n+v+1
                        fn=(-1)**(n+v+ids)
                        isn=in0+ms
                        isf=in0-ms
                        iwv=iv0+mw
                        iwf=iv0-mw
                        nvs=n+v-max(iabs(n-v),iabs(ids))+1
                        do ii=1,nvs
                           it=it-1
                           do ip=1,2
                              do iq=1,2
                                 A=pct(i,j,iq,ip,iwf,isf)
                                 A=fn*wsdt(n,v,ms,it)*A
                                 fhas(it,ip,iq)=
     +                              fhas(it,ip,iq)+A
                                 B=pct(j,i,ip,iq,isn,iwv)
                                 B=wsdt(n,v,ms,it)*B
                                 fnhs(it,ip,iq)=
     +                              fnhs(it,ip,iq)+B
                              enddo
                           enddo
                        enddo
 1901                continue
                     itmin=max(iabs(ids),iabs(n-v))
                     do it=itmin,n+v
                        do m=0,n
                           if(iabs(m-1).gt.it) then
                              cwmf1=0.d0
                           else
                              cwmf1=wmf1(n,v,m,it)
                           endif
                           if(m+1.gt.it) then
                              cwm1=0.d0
                           else
                              cwm1=wm1(n,v,m,it)
                           endif
                           do ip=1,2
                              do iq=1,2
                                 cq=(-1)**iq
                                 A=cwmf1*fhas(it,ip,iq)
                                 A=dconjg(A)
                                 B=cv*A
                                 fhmf1(it,n,m,ip)=
     +                              fhmf1(it,n,m,ip)+A
                                 fhmf1vq(it,n,m,ip)=
     +                              fhmf1vq(it,n,m,ip)+cq*B
                                 A=cwm1*fhas(it,ip,iq)
                                 A=dconjg(A)
                                 B=cv*A
                                 fhm1v(it,n,m,ip)=
     +                              fhm1v(it,n,m,ip)+B
                                 fhm1q(it,n,m,ip)=
     +                              fhm1q(it,n,m,ip)+cq*A
                                 A=cwmf1*fnhs(it,ip,iq)
                                 B=cv*A
                                 fmf1(it,n,m,ip)=
     +                              fmf1(it,n,m,ip)+A
                                 fmf1vq(it,n,m,ip)=
     +                              fmf1vq(it,n,m,ip)+cq*B
                                 A=cwm1*fnhs(it,ip,iq)
                                 B=cv*A
                                 fm1v(it,n,m,ip)=
     +                              fm1v(it,n,m,ip)+B
                                 fm1q(it,n,m,ip)=
     +                              fm1q(it,n,m,ip)+cq*A
                              enddo
                           enddo
                        enddo
                     enddo
19004             continue
19003          continue
               do 1904 n=1,nmax0
                  if(n.gt.nmax(i)) goto 1904
                  itmin=iabs(ids)
                  if(itmin.gt.n+nmax0) goto 1904
                  do 1905 jn=1,nmax0
                     if(jn.gt.nmax(j)) goto 1905
                     if(itmin.gt.jn+nmax0) goto 1905
                     itmax=min(n,jn)+nmax0
                     do 1906 it=itmin,itmax
                        gt=2*it+1
                        do m=0,min(n,jn)
                           do ip=1,2
                              do jp=1,2
                                 A=fhmf1(it,n,m,ip)
                                 B=fmf1(it,jn,m,jp)
                                 A0=A*B
                                 A=fhm1q(it,n,m,ip)
                                 B=fm1q(it,jn,m,jp)
                                 A=gt*(A0+A*B)
                                 A1m(m,n,jn,ip,jp)=
     +                              A1m(m,n,jn,ip,jp)+A
                                 fn=(-1)**(ip+jp+n+jn)
                                 A=fhmf1vq(it,n,m,ip)
                                 B=fmf1vq(it,jn,m,jp)
                                 A0=A*B
                                 A=fhm1v(it,n,m,ip)
                                 B=fm1v(it,jn,m,jp)
                                 A=fn*gt*(A0+A*B)
                                 A2m(m,n,jn,ip,jp)=
     +                              A2m(m,n,jn,ip,jp)+A
                                 if(m.gt.0) goto 19061
                                 temp=((-1)**ip)*gt
                                 A=fhm1q(it,n,0,ip)
                                 B=fmf1(it,jn,2,jp)
                                 A=temp*A*B
                                 B11n(n,jn,ip,jp)=
     +                              B11n(n,jn,ip,jp)+A
                                 nj=ip+n+it
                                 temp=((-1)**nj)*gt
                                 A=fhmf1vq(it,n,1,ip)
                                 B=fmf1(it,jn,1,jp)
                                 A=temp*A*B
                                 B12n(n,jn,ip,jp)=
     +                              B12n(n,jn,ip,jp)+A
                                 nj=ip+n+jn
                                 temp=((-1)**nj)*gt
                                 A=fhmf1vq(it,n,2,ip)
                                 B=fm1v(it,jn,0,jp)
                                 A=temp*A*B
                                 B13n(n,jn,ip,jp)=
     +                              B13n(n,jn,ip,jp)+A
                                 nj=jp+n+jn
                                 temp=((-1)**nj)*gt
                                 A=fhm1v(it,n,0,ip)
                                 B=fmf1vq(it,jn,2,jp)
                                 A=temp*A*B
                                 B21n(n,jn,ip,jp)=
     +                              B21n(n,jn,ip,jp)+A
                                 nj=jp+jn+it
                                 temp=((-1)**nj)*gt
                                 A=fhmf1(it,n,1,ip)
                                 B=fmf1vq(it,jn,1,jp)
                                 A=temp*A*B
                                 B22n(n,jn,ip,jp)=
     +                              B22n(n,jn,ip,jp)+A
                                 temp=((-1)**jp)*gt
                                 A=fhmf1(it,n,2,ip)
                                 B=fm1q(it,jn,0,jp)
                                 A=temp*A*B
                                 B23n(n,jn,ip,jp)=
     +                              B23n(n,jn,ip,jp)+A
19061                            continue
                              enddo
                           enddo
                        enddo
 1906                continue
 1905             continue
 1904          continue
            enddo
         enddo
 1900 continue
      do 1910 ia=1,nang
         iang=2*nang-ia
         dang(ia)=sang*dble(ia-1)
         dang(iang)=180.0d0-dang(ia)                  
         theta=dang(ia)*pione/180.0d0
         xt=dcos(theta)
         call tipitaud(nmax0,xt)
         do i=1,4
            do j=1,4
               mue(i,j,ia)=0.d0
               mue(i,j,iang)=0.d0
            enddo
         enddo
         do i=1,2
            do j=1,2
               do ik=1,2
                  do jk=1,2
                     A0p(i,j,ik,jk)=ci0
                     A1p(i,j,ik,jk)=ci0
                     A0pg(i,j,ik,jk)=ci0
                     A1pg(i,j,ik,jk)=ci0
                     B0p(i,j,ik,jk)=ci0
                     B1p(i,j,ik,jk)=ci0
                     B0pg(i,j,ik,jk)=ci0
                     B1pg(i,j,ik,jk)=ci0
                  enddo
               enddo
            enddo
         enddo
         do n=1,nmax0
            gn=(-1)**n
            itau0=(n-1)*(n+2)/2+1
            tau0p(1)=0.5d0*tau(itau0)
            tau0p(2)=0.d0
            tau0pg(1)=-gn*tau0p(1)
            tau0pg(2)=0.d0
            itau=itau0+1
            gn=-gn
            tau1p(1)=tau(itau)
            tau1p(2)=pi(itau)
            tau1pg(1)=-gn*tau1p(1)
            tau1pg(2)=gn*tau1p(2)
            itau=itau0+2
            gn=-gn
            tau2p(1)=tau(itau)
            tau2p(2)=pi(itau)
            tau2pg(1)=-gn*tau2p(1)
            tau2pg(2)=gn*tau2p(2)
            do jn=1,nmax0
               jtau0=(jn-1)*(jn+2)/2+1
               gn=(-1)**jn
               tau0pj(1)=0.5d0*tau(jtau0)
               tau0pj(2)=0.d0
               tau0pjg(1)=-gn*tau0pj(1)
               tau0pjg(2)=0.d0
               jtau=jtau0+1
               gn=-gn
               tau1pj(1)=tau(jtau)
               tau1pj(2)=pi(jtau)
               tau1pjg(1)=-gn*tau1pj(1)
               tau1pjg(2)=gn*tau1pj(2)
               jtau=jtau0+2
               gn=-gn
               tau2pj(1)=tau(jtau)
               tau2pj(2)=pi(jtau)
               tau2pjg(1)=-gn*tau2pj(1)
               tau2pjg(2)=gn*tau2pj(2)
               do ik=1,2
                  do jk=1,2
                     tau20(ik,jk)=tau2p(ik)*tau0pj(jk)
                     tau11(ik,jk)=tau1p(ik)*tau1pj(jk)
                     tau02(ik,jk)=tau0p(ik)*tau2pj(jk)
                     tau20g(ik,jk)=tau2pg(ik)*tau0pjg(jk)
                     tau11g(ik,jk)=tau1pg(ik)*tau1pjg(jk)
                     tau02g(ik,jk)=tau0pg(ik)*tau2pjg(jk)
                  enddo
               enddo
               do ip=1,2
                  do jp=1,2
                     A=B11n(n,jn,ip,jp)
                     B=B12n(n,jn,ip,jp)
                     cmz=B13n(n,jn,ip,jp)
                     Aj=B21n(n,jn,ip,jp)
                     Bj=B22n(n,jn,ip,jp)
                     cmzj=B23n(n,jn,ip,jp)
                     do ik=1,2
                        do jk=1,2
                           A1=A*tau02(ik,jk)
                           A2=B*tau11(ik,jk)
                           Aj2=cmz*tau20(ik,jk)
                           B1=A1-A2+Aj2                           
                           A1=Aj*tau02(ik,jk)
                           A2=Bj*tau11(ik,jk)
                           Aj2=cmzj*tau20(ik,jk)
                           B2=A1-A2+Aj2
                           A0=B1+B2
                           B0p(ip,jp,ik,jk)=B0p(ip,jp,ik,jk)+A0
                           A0=B1-B2
                           B1p(ip,jp,ik,jk)=B1p(ip,jp,ik,jk)+A0
                           A1=A*tau02g(ik,jk)
                           A2=B*tau11g(ik,jk)
                           Aj2=cmz*tau20g(ik,jk)
                           B1=A1-A2+Aj2
                           A1=Aj*tau02g(ik,jk)
                           A2=Bj*tau11g(ik,jk)
                           Aj2=cmzj*tau20g(ik,jk)
                           B2=A1-A2+Aj2
                           A0=B1+B2
                           B0pg(ip,jp,ik,jk)=B0pg(ip,jp,ik,jk)+A0
                           A0=B1-B2
                           B1pg(ip,jp,ik,jk)=B1pg(ip,jp,ik,jk)+A0
                        enddo
                     enddo
                  enddo
               enddo
               do m=0,min(n,jn)
                  itau=itau0+m
                  gmn=(-1)**(m+n)
                  taup(1)=tau(itau)
                  taup(2)=pi(itau)
                  if(m.eq.0) then
                     taup(1)=0.5d0*taup(1)
                     taup(2)=0.d0
                  endif
                  taupg(1)=-gmn*taup(1)
                  taupg(2)=gmn*taup(2)
                  gmnj=(-1)**(jn+m)
                  jtau=jtau0+m
                  taupj(1)=tau(jtau)
                  taupj(2)=pi(jtau)
                  if(m.eq.0) then
                     taupj(1)=0.5d0*taupj(1)
                     taupj(2)=0.d0
                  endif
                  taupjg(1)=-gmnj*taupj(1)
                  taupjg(2)=gmnj*taupj(2)
                  do ik=1,2
                     do jk=1,2
                        taum(ik,jk)=taup(ik)*taupj(jk)
                        taumg(ik,jk)=taupg(ik)*taupjg(jk)
                     enddo
                  enddo            
                  do ip=1,2
                     do jp=1,2
                        A=A1m(m,n,jn,ip,jp)
                        Aj=A
                        B=A2m(m,n,jn,ip,jp)
                        A=A+B
                        Aj=Aj-B
                        do ik=1,2
                           do jk=1,2
                              B0=A*taum(ik,jk)
                              B1=Aj*taum(ik,jk)
                              A0p(ip,jp,ik,jk)=
     +                           A0p(ip,jp,ik,jk)+B0
                              A1p(ip,jp,ik,jk)=
     +                           A1p(ip,jp,ik,jk)+B1
                              B0=A*taumg(ik,jk)
                              B1=Aj*taumg(ik,jk)
                              A0pg(ip,jp,ik,jk)=
     +                           A0pg(ip,jp,ik,jk)+B0
                              A1pg(ip,jp,ik,jk)=
     +                           A1pg(ip,jp,ik,jk)+B1
                           enddo
                        enddo
                     enddo
                  enddo
               enddo
            enddo
         enddo         
         do ip=1,2
            do jp=1,2
               temp=A0p(ip,jp,ip,jp)+A0p(ip,jp,3-ip,3-jp)
               mue(1,1,ia)=mue(1,1,ia)+temp
               temp=B0p(ip,jp,3-ip,3-jp)-B0p(ip,jp,ip,jp)
               mue(1,2,ia)=mue(1,2,ia)+temp
               A0=A1p(ip,jp,ip,jp)+A1p(ip,jp,3-ip,3-jp)
               A0=A0-B1p(ip,jp,ip,jp)+B1p(ip,jp,3-ip,3-jp)
               temp=-dimag(A0)
               mue(1,3,ia)=mue(1,3,ia)+temp
               temp=-A0
               mue(1,4,ia)=mue(1,4,ia)+temp
               temp=A0p(ip,jp,ip,jp)-A0p(ip,jp,3-ip,3-jp)
               mue(2,1,ia)=mue(2,1,ia)+temp
               temp=B0p(ip,jp,ip,jp)+B0p(ip,jp,3-ip,3-jp)
               mue(2,2,ia)=mue(2,2,ia)-temp
               A0=A1p(ip,jp,ip,jp)+A1p(ip,jp,3-ip,3-jp)
               A0=A0+B1p(ip,jp,ip,jp)-B1p(ip,jp,3-ip,3-jp)
               temp=-dimag(A0)
               mue(2,3,ia)=mue(2,3,ia)+temp
               temp=-A0
               mue(2,4,ia)=mue(2,4,ia)+temp
               temp=0.5d0*dimag(A1p(ip,jp,ip,3-jp))
               mue(3,1,ia)=mue(3,1,ia)+temp
               temp=-0.5d0*dimag(B1p(ip,jp,ip,3-jp))
               mue(3,2,ia)=mue(3,2,ia)+temp
               A0=A0p(ip,jp,ip,3-jp)-A0p(ip,jp,3-ip,jp)
               A0=A0-B0p(ip,jp,ip,3-jp)-B0p(ip,jp,3-ip,jp)
               temp=A0
               mue(3,3,ia)=mue(3,3,ia)+temp
               temp=-dimag(A0)
               mue(3,4,ia)=mue(3,4,ia)+temp
               temp=-0.5d0*A1p(ip,jp,ip,3-jp)
               mue(4,1,ia)=mue(4,1,ia)+temp
               temp=0.5d0*B1p(ip,jp,ip,3-jp)
               mue(4,2,ia)=mue(4,2,ia)+temp
               A0=A0p(ip,jp,ip,3-jp)+A0p(ip,jp,3-ip,jp)
               A0=A0-B0p(ip,jp,ip,3-jp)+B0p(ip,jp,3-ip,jp)
               temp=dimag(A0)
               mue(4,3,ia)=mue(4,3,ia)+temp
               temp=A0
               mue(4,4,ia)=mue(4,4,ia)+temp
               if(ia.eq.iang) goto 19101
               temp=A0pg(ip,jp,ip,jp)+A0pg(ip,jp,3-ip,3-jp)
               mue(1,1,iang)=mue(1,1,iang)+temp
               temp=B0pg(ip,jp,3-ip,3-jp)-B0pg(ip,jp,ip,jp)
               mue(1,2,iang)=mue(1,2,iang)+temp
               A0=A1pg(ip,jp,ip,jp)+A1pg(ip,jp,3-ip,3-jp)
               A0=A0-B1pg(ip,jp,ip,jp)+B1pg(ip,jp,3-ip,3-jp)
               temp=-dimag(A0)
               mue(1,3,iang)=mue(1,3,iang)+temp
               temp=-A0
               mue(1,4,iang)=mue(1,4,iang)+temp
               temp=A0pg(ip,jp,ip,jp)-A0pg(ip,jp,3-ip,3-jp)
               mue(2,1,iang)=mue(2,1,iang)+temp
               temp=B0pg(ip,jp,ip,jp)+B0pg(ip,jp,3-ip,3-jp)
               mue(2,2,iang)=mue(2,2,iang)-temp
               A0=A1pg(ip,jp,ip,jp)+A1pg(ip,jp,3-ip,3-jp)
               A0=A0+B1pg(ip,jp,ip,jp)-B1pg(ip,jp,3-ip,3-jp)
               temp=-dimag(A0)
               mue(2,3,iang)=mue(2,3,iang)+temp
               temp=-A0
               mue(2,4,iang)=mue(2,4,iang)+temp
               temp=0.5d0*dimag(A1pg(ip,jp,ip,3-jp))
               mue(3,1,iang)=mue(3,1,iang)+temp
               temp=-0.5d0*dimag(B1pg(ip,jp,ip,3-jp))
               mue(3,2,iang)=mue(3,2,iang)+temp
               A0=A0pg(ip,jp,ip,3-jp)-A0pg(ip,jp,3-ip,jp)
               A0=A0-B0pg(ip,jp,ip,3-jp)-B0pg(ip,jp,3-ip,jp)
               temp=A0
               mue(3,3,iang)=mue(3,3,iang)+temp
               temp=-dimag(A0)
               mue(3,4,iang)=mue(3,4,iang)+temp
               temp=-0.5d0*A1pg(ip,jp,ip,3-jp)
               mue(4,1,iang)=mue(4,1,iang)+temp
               temp=0.5d0*B1pg(ip,jp,ip,3-jp)
               mue(4,2,iang)=mue(4,2,iang)+temp
               A0=A0pg(ip,jp,ip,3-jp)+A0pg(ip,jp,3-ip,jp)
               A0=A0-B0pg(ip,jp,ip,3-jp)+B0pg(ip,jp,3-ip,jp)
               temp=dimag(A0)
               mue(4,3,iang)=mue(4,3,iang)+temp
               temp=A0
               mue(4,4,iang)=mue(4,4,iang)+temp
19101          continue
            enddo
         enddo
         temp=mue(1,1,ia)+mue(1,2,ia)+mue(2,1,ia)+mue(2,2,ia)
         i22(ia)=0.5d0*temp
         temp=mue(1,1,ia)-mue(1,2,ia)-mue(2,1,ia)+mue(2,2,ia)
         i11(ia)=0.5d0*temp
         temp=mue(1,1,ia)+mue(1,2,ia)-mue(2,1,ia)-mue(2,2,ia)
         i21(ia)=0.5d0*temp
         temp=mue(1,1,ia)-mue(1,2,ia)+mue(2,1,ia)-mue(2,2,ia)
         i12(ia)=0.5d0*temp
         if(ia.eq.iang) goto 1910
         temp=mue(1,1,iang)+mue(1,2,iang)+mue(2,1,iang)+mue(2,2,iang)
         i22(iang)=0.5d0*temp
         temp=mue(1,1,iang)-mue(1,2,iang)-mue(2,1,iang)+mue(2,2,iang)
         i11(iang)=0.5d0*temp
         temp=mue(1,1,iang)+mue(1,2,iang)-mue(2,1,iang)-mue(2,2,iang)
         i21(iang)=0.5d0*temp
         temp=mue(1,1,iang)-mue(1,2,iang)+mue(2,1,iang)-mue(2,2,iang)
         i12(iang)=0.5d0*temp
 1910 continue
      cbak=i11(2*nang-1)
      do i=1,nang2
         inat(i)=i11(i)+i22(i)+i12(i)+i21(i)
         pol(i)=(i11(i)+i12(i)-i22(i)-i21(i))/inat(i)
      enddo
 2000 cz=twopi/k**2
      csca=csca*cz
      cext=cext*cz
      cabs=cext-csca
      cpr=cpr*cz
      assym=cpr/csca
      cbak=2.d0*cbak*cz
      write(6,'(5x,a6,7x,a6,7x,a6,7x,a6,7x,a5,6x,a12)')
     +'<Cext>','<Cabs>','<Csca>','<Cbak>','<Cpr>','<cos(theta)>'
      write(6,'(2x,6e13.5)') cext,cabs,csca,cbak,cext-cpr,assym
      cscax=0.d0
      cextx=0.d0
      cprx=0.d0
      do i=1,nL
         cscax=cscax+cscai(i)
         cextx=cextx+cexti(i)
         cprx=cprx+cpri(i)
      enddo
      cscax=cscax*cz
      cextx=cextx*cz
      cprx=cprx*cz
      cabsx=cextx-cscax
      assymx=cprx/cscax
      cbakx=cbak	
      write(6,'(2x,6e13.5)') cextx,cabsx,cscax,cbakx,cextx-cprx,assymx
      do i=1,nL
         cexti(i)=cexti(i)*cz
         cscai(i)=cscai(i)*cz
         cabsi(i)=cexti(i)-cscai(i)
         cpri(i)=cpri(i)*cz
         assymi(i)=cpri(i)/csca
         write(6,'(i5,5e15.6)') i,cexti(i),cabsi(i),
     +                    cscai(i),cpri(i),assymi(i)
      enddo
c      flout='cr'//fileout
c      open(33,file=flout,status='unknown')
c      write(33,'(a20,a47)') 
c     +     flout,'(Total and individual-particle cross sections)'
c      write(33,'(a32,2x,a22)') 
c     +     'input sphere-aggregate filename:',FLNAME
c      write(33,'(12x,a4,11x,a4,11x,a4,11x,a3,8x,a12)')
c     +    'Cext','Cabs','Csca','Cpr','<cos(theta)>'
c      write(33,'(a5,5e15.6)') 'total',cext,cabs,csca,cext-cpr,
c     +     assym
c      do i=1,nL
c         write(33,'(i5,5e15.6)') i,cexti(i),cabsi(i),cscai(i),
c     +           cpri(i),assymi(i)
c      enddo
c      close(33)       
      cz=pione*gcvr*gcvr
      cscav=csca/cz
      cextv=cext/cz
      cprv=cpr/cz
      cbakv=cbak/cz
      cabsv=cextv-cscav
      temp=gcvr*gcvr/gcs
      cscas=cscav*temp
      cexts=cextv*temp
      cabss=cabsv*temp
      cbaks=cbakv*temp
      cprs=cprv*temp
 222  format(6e13.5)
 221  format(3x,a7,6x,a7,6x,a7,6x,a7,6x,a6,5x,a12)        
      write(6,221)
     +'<Qextv>','<Qabsv>','<Qscav>','<Qbakv>','<Qprv>','<cos(theta)>'	
      write(6,222) cextv,cabsv,cscav,cbakv,cprv,assym
      write(6,221)
     +'<Qexts>','<Qabss>','<Qscas>','<Qbaks>','<Qprs>','<cos(theta)>'
      write(6,222) cexts,cabss,cscas,cbaks,cprs,assym
      open(12,file=fileout,status='unknown')
      write(12,'(/)')
      write(12,'(1x,a15,a13,a18,a4,f8.3,a5,f8.3)') fileout, 
     +        ' input file: ',FLNAME,' xv:',xv,'  xs:',xs
      if(idscmt.lt.0) then
         write(12,'(/)')
         write(12,*) '*** backscattering and scatttering matrix',
     +                 ' are not calculated ***'
      endif
      write(12,'(/)')
      write(12,221)
     +'<Cext>','<Cabs>','<Csca>','<Cbak>','<Cpr>','<cos(theta)>'
      write(12,222) 
     +   cext,cabs,csca,cbak,cext-cpr,assym
      write(12,221)
     +'<Qextv>','<Qabsv>','<Qscav>','<Qbakv>','<Qprv>','<cos(theta)>'
      write(12,222) cextv,cabsv,cscav,cbakv,cprv,assym
      write(12,221)
     +'<Qexts>','<Qabss>','<Qscas>','<Qbaks>','<Qprs>','<cos(theta)>'
      write(12,222) cexts,cabss,cscas,cbaks,cprs,assym
      if(idscmt.lt.0) goto 2001
      write(12,'(/)')
      write(12,'(2x,a4,4x,a7,4x,a6,4x,a7,6x,a7,6x,a7,6x,a7)')
     +'s.a.','<total>','<pol.>','<S1*S1>','<S4*S4>','<S3*S3>','<S2*S2>'               
      do i=1,nang2
         write(12,'(f6.1,e13.5,f8.4,4e13.5)') 
     +      dang(i),inat(i),pol(i),i11(i),i21(i),i12(i),i22(i) 
      enddo
      write(12,'(/)')	
      write(12,'(1x,a50)') 
     +   'Scattering matrix (4X4 for each scattering angle):'
      do i=1,nang2
         write(12,'(f7.1,4e16.7)') 
     +      dang(i),mue(1,1,i),mue(1,2,i),mue(1,3,i),mue(1,4,i)
         write(12,'(7x,4e16.7)') 
     +              mue(2,1,i),mue(2,2,i),mue(2,3,i),mue(2,4,i)
         write(12,'(7x,4e16.7)') 
     +              mue(3,1,i),mue(3,2,i),mue(3,3,i),mue(3,4,i)
         write(12,'(7x,4e16.7)') 
     +              mue(4,1,i),mue(4,2,i),mue(4,3,i),mue(4,4,i)
      enddo
 2001 close(12)
      STOP
      END

      SUBROUTINE abMiexud(X,REFREL,NP,NMAX,NM,AN,BN,NADD,
     +                    RSR,RSI,RSX,PX,AR,AI,BR,BI,EPS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      integer NMAX,NM,NADD,NSTOP,NX,K,N
      COMPLEX*16 REFREL,AN(NP),BN(NP)
      DOUBLE PRECISION AR(NP),AI(NP),BR(NP),BI(NP),
     +   RSR(NMAX),RSI(NMAX),RSX(NMAX),PX(NMAX)
      IF(EPS.GT.1.D0.OR.EPS.LT.0.D0) EPS=1.0D-20
      IF(NADD.NE.0) EPS=0.D0
      CTC=EPS
      XM=DBLE(REFREL)
      YM=DIMAG(REFREL)
      XMX=X*XM
      YMX=X*YM
      RP2=XMX*XMX+YMX*YMX
      NSTOP=X+4.D0*X**.3333D0
      NSTOP=NSTOP+2
      NM=NSTOP+NADD
      XN=DSQRT(XM**2+YM**2)*X
      NX=1.1D0*XN+10.D0
      if(NX-NM.lt.10) NX=NM+10
      write(6,*) 'Wiscombe criterion:',NSTOP 
      write(6,*) 'NADD:',NADD
      write(6,*) 'NX:',NX
      IF(NX.GT.NMAX) THEN
         WRITE(6,*) 'Parameter NXMAX too small'
         WRITE(6,*) '  NXMAX must be greater than', NX
         WRITE(6,*) 'Please correct NXMAX in main code,' 
         WRITE(6,*) '  recompile, then try again'
         STOP
      ENDIF
      IF(NM.GT.NP) THEN
         WRITE(6,*) 'Parameter np too small'
         WRITE(6,*) '  np must be greater than', NM
         WRITE(6,*) 'Please correct np in gmm01f.par,' 
         WRITE(6,*) '  recompile the code, then try again'
         STOP
      ENDIF
C  DOWN RECURSION FOR RATIOS RSR,RSI,RSX,PNR,PNI,PX
      PNX=X/DBLE(2*NX+3)
      PNR=XMX/DBLE(2*NX+3)
      PNI=YMX/DBLE(2*NX+3)
      DO 5 K=1,NX
         N=NX-K+1
         CN=DBLE(N)
         ALN=(2.D0*CN+1.D0)*XMX/RP2-PNR
         BEN=(2.D0*CN+1.D0)*YMX/RP2+PNI
         RSR(N)=-CN*XMX/RP2+ALN
         RSI(N)=CN*YMX/RP2-BEN
         PZD=ALN*ALN+BEN*BEN
         PNR=ALN/PZD
         PNI=BEN/PZD
         RSX(N)=(CN+1.D0)/X-PNX
         IF(N.EQ.1) GO TO 20
         PNX=X/(2.D0*CN+1.D0-PNX*X)
         PX(N)=PNX
    5    CONTINUE
   20 SNM1X=DSIN(X)
      CNM1X=DCOS(X)
      IF(X-0.1D0)21,22,22
   21 SNX=X**2./3.D0-X**4./30.D0+X**6./840.D0-X**8./45360.D0
      GO TO 23
   22 SNX=SNM1X/X-CNM1X
   23 CNX=CNM1X/X+SNM1X
      DO 10 N=1,NX
         PX(N)=SNX      !preparing for the calculation of Cabs
         C=DBLE(N)
         DCNX=CNM1X-C*CNX/X
         DSNX=RSX(N)*SNX
C  CALCULATION OF EXTERIOR COEFFICIENTS AN AND BN
         ANNR=RSR(N)*SNX-XM*DSNX
         ANNI=RSI(N)*SNX-YM*DSNX
         TA1=RSR(N)*SNX-RSI(N)*CNX
         TA2=RSI(N)*SNX+RSR(N)*CNX
         ANDR=TA1-XM*DSNX+YM*DCNX
         ANDI=TA2-XM*DCNX-YM*DSNX
         AND=ANDR*ANDR+ANDI*ANDI
         BNNR=(XM*RSR(N)-YM*RSI(N))*SNX-DSNX
         BNNI=(XM*RSI(N)+YM*RSR(N))*SNX
         TB1=RSR(N)*SNX-RSI(N)*CNX
         TB2=RSR(N)*CNX+RSI(N)*SNX
         BNDR=XM*TB1-YM*TB2-DSNX
         BNDI=XM*TB2+YM*TB1-DCNX
         BND=BNDR*BNDR+BNDI*BNDI
         AR(N)=(ANNR*ANDR+ANNI*ANDI)/AND
         AI(N)=(ANNI*ANDR-ANNR*ANDI)/AND
         BR(N)=(BNNR*BNDR+BNNI*BNDI)/BND
         BI(N)=(BNNI*BNDR-BNNR*BNDI)/BND
C  MIE SERIES CONVERGENCE TEST IS MADE BY TESTING AN'S AND BN'S
         TI=AR(N)*AR(N)+AI(N)*AI(N)+BR(N)*BR(N)+BI(N)*BI(N)
         TI=TI/(AR(1)*AR(1)+AI(1)*AI(1)+BR(1)*BR(1)+BI(1)*BI(1))
         IF(TI-CTC) 16,18,18
   18    IF(NM-N) 15,15,6
    6    IF(N-NX)7,15,15
    7    M=N+1
         SNX=PX(M)*SNX
         CNM2X=CNM1X
         CNM1X=CNX
         CNX=(2.D0*C+1.D0)*CNM1X/X-CNM2X
   10 CONTINUE
      GO TO 15
   16 WRITE(6,*) '*** NOTE THAT THE FIELD-EXPANSION TRANCATION'
      WRITE(6,*) '*** IS DETERMINED BY eps GIVEN IN THE INPUT'
      WRITE(6,*) '*** FILE gmm01f.in'
      WRITE(6,*) '*** IN CASE YOU NEED A HIGHER ORDER, eps MUST'
      WRITE(6,'(a,e9.1)')
     +          ' *** BE SMALLER THAN THE CURRENT VALUE',EPS
   15 NM=N
      DO I=1,NM
         AN(I)=DCMPLX(AR(I),-AI(I))
         BN(I)=DCMPLX(BR(I),-BI(I))
      ENDDO
      RETURN
      END

	subroutine cofsrd(nmax)
	include 'gmm01f.par'
	parameter (nmp=np*(np+2))	
	double precision cofsr(nmp),lnfacd,c
	common/crot/cofsr
        i=0
	do n=1,nmax
	   do m=-n,n
	      i=i+1       
	      c=lnfacd(dble(n-m))-lnfacd(dble(n+m))
	      cofsr(i)=0.5d0*c
c	      c=0.5d0*c
c	      cofsr(i)=dexp(c)
           enddo
        enddo
	return
	end

	subroutine cofd0(nmax)
	implicit double precision (a-h,o-z)	
	include 'gmm01f.par'
	parameter (nmp=np*(np+2))	
	parameter (ni0=np*(np+1)*(2*np+1)/3+np*np)	
	integer v
	double precision lnfacd
	common/cofmnv0/cof0(ni0)
	common/crot/cofsr(nmp)	
        common/fnr/fnr(0:4*(np+1))
	i=0
	sm=-0.5d0*dble((-1)**nmax)
	do m=-nmax,nmax
	   ns=max(1,iabs(m))
	   sm=-sm
	   do n=ns,nmax
	      inm=n*(n+1)-m
	      do v=ns,nmax
	         i=i+1
	         ivm=v*(v+1)+m
	         c=cofsr(inm)+cofsr(ivm)	            
	         c=sm*dexp(c)
	         c0=fnr(2*n+1)*fnr(2*v+1)
	         c1=fnr(n)*fnr(v)*fnr(n+1)*fnr(v+1)
	         c0=c0/c1
                 cof0(i)=c*c0
              enddo
           enddo
        enddo
	return
	end

	subroutine cofnv0(nmax)
	include 'gmm01f.par'
	integer n,v
	double precision c1,lnfacd,cnv(np,np)
	common/cfnv/cnv
	do n=1,nmax
	   do v=n,nmax
	      c1=lnfacd(dble(2*n))+lnfacd(dble(2*v))
	      c1=c1-lnfacd(dble(2*n+2*v))
	      c1=c1+2.d0*lnfacd(dble(n+v))
	      c1=c1-lnfacd(dble(n))-lnfacd(dble(v))
	      cnv(n,v)=c1
	   enddo
	enddo
	return
	end

C  subroutine gau0.f generates tabulated values for  
C  Gaunt coefficients up to n=v=n_max
	subroutine gau0(nmax)
	include 'gmm01f.par'
	parameter (ni0=np*(np+1)*(2*np+1)/3+np*np)
	parameter (ng0=np*(2*np**3+10*np**2+19*np+5)/6)
	integer v,qmax,uvmax,iga0(ni0)
	double precision ga0(ng0)
	common/g0/ga0
	common/ig0/iga0
        na=0
        uvmax=nmax*(nmax+2)
        i=0
	do m=-nmax,nmax
	   ns=max(1,iabs(m))
	   do n=ns,nmax
	      do v=ns,nmax
	         call gxurcd0(-m,n,v,qmax,na)
	         i=i+1
                 iga0(i)=na
	         na=na+qmax+1
              enddo
	   enddo
        enddo
	return
	end

c  transforms the rectangular coordinates (x,y,z)
c  to spherical coordinates (r,theta,phi)
	subroutine carsphd(x,y,z,r,xt,sphi,cphi)
	double precision x,y,z,r,xt,sphi,cphi
	r=dsqrt(x*x+y*y+z*z)
	if(r.eq.0.d0) then
	   xt=1.d0
	   sphi=0.d0
	   cphi=1.d0
	   return
	endif
	xt=z/r
	if(y.eq.0.d0.and.x.eq.0.d0) then
	   sphi=0.d0
	   cphi=1.d0
	   return
	endif
	sphi=dsqrt(x*x+y*y)
	cphi=x/sphi
	sphi=y/sphi
	return
	end

c  subroutine besseljd.f  (in double precision arithmetic)
c  returns an array of the spherical Bessel function of the  
c  first kind with a real argument: j_0,j_1,j_2,...,j_{NC}
c  uses Ru Wang's ratio method for the downward recursive 
c  calculation of the Riccati-Bessel function Psi_n(z)=z j_n(z) 
c  [see Xu et al., Physical Review E, v.60, 2347-2365 (1999)]
      SUBROUTINE besseljd(NC,X,BESJ)
      INTEGER NC,NX,K,N
      DOUBLE PRECISION X,BESJ(0:NC),PN,CN,X2
      DO K=1,NC
         BESJ(K)=0.D0
      ENDDO
      IF(DABS(X).LT.1.D-12) THEN
         BESJ(0)=1.D0
         BESJ(1)=1.D0/3.D0*X
         RETURN
      ENDIF
c  down-recursively calculating an array of the ratio functions  
c  P_{NC},P_{NC-1},...,P(2) stored in the same array for j_n, 
c  starting with an asymptotic value P_{NX+1}=X/(2 NX+3) at the 
c  highest order NX+1, where NX=NC+1.1X+10 
      NX=1.1D0*X+10.D0
      NX=NC+NX   
      PN=X/DBLE(2*NX+3)
      DO 5 K=1,NX-1
         N=NX-K+1
         CN=DBLE(N)
         PN=X/(DBLE(2*N+1)-PN*X)
         IF(N.GT.NC) GOTO 5
         BESJ(N)=PN
    5 CONTINUE
C  calculating j_0(x) and j_1(x)
      IF(DABS(X)-0.1D0) 10,11,11
   10 X2=X*X
      BESJ(0)=1.D0-X2/72.D0
      BESJ(0)=1.D0-X2/42.D0*BESJ(0)
      BESJ(0)=1.D0-X2/20.D0*BESJ(0)
      BESJ(0)=1.D0-X2/6.D0*BESJ(0)
      BESJ(1)=1.D0/45360.D0-X2/3991680.D0
      BESJ(1)=1.D0/840.D0-X2*BESJ(1)
      BESJ(1)=1.D0/30.0d0-X2*BESJ(1)
      BESJ(1)=X*(1.D0/3.0d0-X2*BESJ(1))
      GOTO 12
   11 BESJ(0)=DSIN(X)/X
      BESJ(1)=(DSIN(X)/X-DCOS(X))/X      
c  calculating j_2,...,j_{NC} 
   12 DO 20 N=2,NC
         BESJ(N)=BESJ(N)*BESJ(N-1)
   20 CONTINUE
      RETURN
      END

c  sub. besselyd.f  (in double precision arithmetic)
c  returns an array of the spherical Bessel function of
c  the second kind with a real argument: y_0,y_1,...,y_n

       subroutine besselyd(n,x,besy)
       integer i,n
       double precision x,besy(0:n),besyn,x2
       if(x.eq.0.d0) then
         write(6,*) 'bad argument in sub. besselyd'
         stop
       endif
       if(dabs(x)-0.1d0)10,11,11
  10   x2=x*x
       besyn=1.d0-x2/72.d0
       besyn=1.d0-x2/42.d0*besyn
       besyn=1.d0-x2/20.d0*besyn
       besyn=1.d0-x2/6.d0*besyn
       besy(0)=1.d0-x2/56.d0
       besy(0)=1.d0-x2/30.d0*besy(0)
       besy(0)=1.d0-x2/12.d0*besy(0)
       besy(0)=1.d0-0.5d0*x2*besy(0)
       besy(0)=-besy(0)/x
       goto 12
  11   besyn=dsin(x)/x
       besy(0)=-dcos(x)/x
  12   besy(1)=besy(0)/x-besyn
       do i=2,n
         besy(i)=dble(2*i-1)/x*besy(i-1)-besy(i-2)
       enddo
       return
       end

c  subroutine rotcoef is originally written by Mackowski, Fuller, and 
c  Mishchenko (taken from the code scsmfo1b.for released to public by 
c  the authors at 8/1999)
c  the rotational coefficients are required for the implementation of 
c  Mackowski's three-step numerical technique in subroutine rtr.f for 
c  decomposition of translation matrix into rotational and axial 
c  translational parts [see Mackowski, Proc. R. Soc. Lond. A 433, 599 
c  (1991)]
c  Yu-lin Xu   12/2000         

      subroutine rotcoef(cbe,nmax)
      include 'gmm01f.par'
      parameter (nmp=np*(np+2))
      implicit double precision (a-h,o-z)
      double precision dk0(-2*np:2*np),dk01(-2*np:2*np)
      common/rot/bcof(0:np+2),dc(-np:np,0:nmp)
      common/fnr/fnr(0:4*(np+1))
      sbe=dsqrt((1.d0+cbe)*(1.d0-cbe))
      cbe2=.5d0*(1.d0+cbe)
      sbe2=.5d0*(1.d0-cbe)
      in=1
      dk0(0)=1.d0
      sben=1.d0
      dc(0,0)=1.d0
      dk01(0)=0.d0
      do n=1,nmax
         nn1=n*(n+1)
         in=-in
         sben=sben*sbe/2.d0
         dk0(n)=dble(in)*sben*bcof(n)
         dk0(-n)=dble(in)*dk0(n)
         dk01(n)=0.d0
         dk01(-n)=0.d0
         dc(0,nn1+n)=dk0(n)
         dc(0,nn1-n)=dk0(-n)
         do k=-n+1,n-1
            kn=nn1+k
            dkt=dk01(k)
            dk01(k)=dk0(k)
            dk0(k)=(cbe*dble(n+n-1)*dk01(k)-fnr(n-k-1)*fnr(n+k-1)*dkt)
     1             /(fnr(n+k)*fnr(n-k))
            dc(0,kn)=dk0(k)
         enddo
         im=1
         do m=1,n
            im=-im
            fmn=1.d0/fnr(n-m+1)/fnr(n+m)
            m1=m-1
            dkm0=0.d0
            do k=-n,n
               kn=nn1+k
               dkm1=dkm0
               dkm0=dc(m1,kn)
               if(k.eq.n) then
                  dkn1=0.d0
               else
                  dkn1=dc(m1,kn+1)
               endif
               dc(m,kn)=(fnr(n+k)*fnr(n-k+1)*cbe2*dkm1
     1           -fnr(n-k)*fnr(n+k+1)*sbe2*dkn1
     1              -dble(k)*sbe*dc(m1,kn))*fmn
               dc(-m,nn1-k)=dble((-1)**(k)*im)*dc(m,kn)
            enddo
         enddo
      enddo
      return
      end

c  subroutine cofxuds0.f returns the two classes of vector 
c  (axial) translation coefficients for a given combination of 
c  (m,n,m,v) and a given dimensionless translation distance kd 
cu uses subroutine gid0.f 
	subroutine cofxuds0(nmax,m,n,v,sja,sya,A,B,Aj,Bj)
	implicit double precision (a-h,o-z)
	include 'gmm01f.par'
	parameter (ni0=np*(np+1)*(2*np+1)/3+np*np)
	parameter (ng0=np*(2*np**3+10*np**2+19*np+5)/6)
	integer v,p,qmax
	double precision sja(0:n+v+1),sya(0:n+v+1)
	complex*16 A,B,Aj,Bj,signz
	common/ig0/iga0(ni0)
	common/g0/ga0(ng0)
	common/cofmnv0/cof0(ni0)
	fa(m,p)=dble(-2*m*p*(p-1))
	fb(n,v,p)=dble(p*p-(n+v+1)*(n+v+1))
     +            *dble(p*p-(n-v)*(n-v))/dble(4*p*p-1)
	if(iabs(m).gt.n.or.iabs(m).gt.v) then 
	   write(6,*) '|m|>n or v in subroutine cofxuds0.f'
	   stop
	endif
	A=0.d0
	B=0.d0
	Aj=0.d0
	Bj=0.d0
	call gid0(nmax,m,n,v,id)
	c=cof0(id)
	ig=iga0(id)
	nv2=v*(v+1)+n*(n+1)
	signz=dcmplx(0.d0,1.d0)**(n+v)
	p=n+v+2
	qmax=min(n,v)
	do i=1,qmax+1
	   p=p-2
	   cp=dble(nv2-p*(p+1))*ga0(ig+i)
	   sj=sja(p)
	   sy=sya(p)
	   A=A+dcmplx(sj,sy)*signz*cp
	   Aj=Aj+sj*signz*cp
	   signz=-signz
	enddo
	A=A*c
	Aj=Aj*c
	if(m.eq.0) return
 	signz=dcmplx(0.d0,1.d0)**(n+v+1)
 	p=n+v
	do 20 i=1,qmax
	   p=p-2
	   signz=-signz
	   if(i.eq.1) then
              cp=dble(2*p+3)*fa(m,p+3)      
              cp=cp*ga0(ig+1)/dble((p+3)*(p+2))
              goto 21
           endif
           if(i.eq.qmax) then
              if(p.eq.0) goto 22
              nv2=p*(p+1)
              cp=dble(2*p+3)*fa(m,p+1)
              cp=-cp*ga0(ig+i+1)/dble(nv2)
              goto 21
           endif
 22	   c4=fa(m,p+2)
           cp=-dble((p+1)*(p+2))*fb(n,v,p+2)*ga0(ig+i)
	   cp=cp+dble((p+2)*(p+1))*fb(n,v,p+1)*ga0(ig+i+1)
 	   cp=cp*dble(2*p+3)/c4
 21	   sj=sja(p+1)
	   sy=sya(p+1)
 	   B=B+dcmplx(sj,sy)*signz*cp
 	   Bj=Bj+sj*signz*cp
 20	continue
	B=B*c
	Bj=Bj*c
	return
	end


c  subroutine tipitaud.f 
c  calculates pi(m,n) & tau(m,n) up to a specified nmax for all 
c  m=0,1,...n at a given x
c  pi(m,n) & tau(m,n) calculated are normalized by  
c         C_mn=[(2n+1)(n-m)!/n/(n+1)/(n+m)!]^(1/2)
c  Yu-lin Xu    12/2000

	subroutine tipitaud(nmax,x)
	include 'gmm01f.par'
	parameter (nmp0=(np+1)*(np+4)/2)
	implicit double precision (a-h,o-z)
	common/fnr/fnr(0:4*(np+1))
	common/pitau/pi(nmp0),tau(nmp0)	

	nt=(nmax+1)*(nmax+4)/2          ! calculates pi up to nmax+1
	if(nt.gt.nmp0.or.dabs(x).gt.1.d0) then
	   write(6,*) 'dimension or argument wrong in sub. tipitaud'
	   write(6,*) 'argument: ',x
	   stop
	endif
	sx=dsqrt(1.d0-x*x)
	pi(1)=0.d0                     ! pi(0,1)  pi(0,n)=0 when m=0
	pi(2)=dsqrt(.75d0)             ! pi(1,1)
	pi(3)=0.d0                     ! pi(0,2)
	t125=dsqrt(1.25d0) 
	pi(4)=t125*x                      ! pi(1,2)  
	pi(5)=t125*sx                     ! pi(2,2)
	imn=5
	do i=3,nmax+1
	   n=i
	   imn=imn+1
	   pi(imn)=0.d0                ! pi(0,n)=0
	   do 11 j=2,n
	      m=j-1
	      imn=imn+1
	      i1=(n-2)*(n+1)/2+m+1
	      if(m.eq.n-1) then
	         pi(imn)=fnr(n-1)*fnr(2*n+1)/fnr(n+1)*x*pi(i1)
	         goto 11
	      endif
	      i2=(n-3)*n/2+m+1
	      t=fnr(n)*fnr(2*n-3)
	      t=fnr(n-2)*fnr(n-m-1)*fnr(n+m-1)/t
	      pi(imn)=fnr(2*n-1)*x*pi(i1)-t*pi(i2)
	      t=fnr(n+1)*fnr(n-m)*fnr(n+m)
	      t=fnr(n-1)*fnr(2*n+1)/t
	      pi(imn)=t*pi(imn)
 11	   continue
	   imn=imn+1
	   i1=(n-2)*(n+1)/2+n
	   t=fnr(n-1)*fnr(n+1)
	   t=dsqrt(.5d0)*fnr(n)*fnr(2*n+1)/t
	   pi(imn)=t*sx*pi(i1)
	enddo
	tx=x*sx
	tau(1)=-dsqrt(1.5d0)*sx          ! tau(0,1)
	tau(2)=pi(2)*x                   ! tau(1,1)
	tau(3)=-dsqrt(7.5d0)*tx          ! tau(0,2)   
	tau(4)=t125*(2.d0*x*x-1.d0)      ! tau(1,2)
	tau(5)=t125*tx                   ! tau(2,2)
	imn=5
	do i=3,nmax
	   n=i
	   do 21 j=1,n+1
	      m=j-1
	      imn=imn+1
	      if(m.eq.0) then
	         i1=(n-2)*(n+1)/2+1
	         i2=(n-3)*n/2+1
	         t=fnr(2*n-3)
	         t=fnr(n-2)*fnr(n)/t
	         tau(imn)=fnr(2*n-1)*x*tau(i1)-t*tau(i2)
	         t=fnr(n-1)*fnr(n+1)
	         t=fnr(2*n+1)/t
	         tau(imn)=t*tau(imn)
	         goto 21
	      endif
	      i1=n*(n+3)/2+m+1
	      t=fnr(n)*fnr(2*n+3)
	      t=fnr(n+2)*fnr(2*n+1)*fnr(n-m+1)*fnr(n+m+1)/t
	      tau(imn)=t*pi(i1)-dble(n+1)*x*pi(imn)
	      tau(imn)=tau(imn)/dble(m)
 21	   continue
	enddo
	return
	end

c  lnfacd.f  (double precision function)
c  returns ln(z!)  z>-1.0
c  based on Lanczos' method [see Xu, Journal of Computational 
c  Physics, v.139, 137-165 (1998)]
      double precision function lnfacd(z)
      integer i
      double precision z,a,b,cp,c0(11)
      data c0/0.16427423239836267d5, -0.48589401600331902d5,
     +        0.55557391003815523d5, -0.30964901015912058d5,
     +        0.87287202992571788d4, -0.11714474574532352d4,
     +        0.63103078123601037d2, -0.93060589791758878d0,
     +        0.13919002438227877d-2,-0.45006835613027859d-8,
     +        0.13069587914063262d-9/ 
      a=1.d0
      cp=2.5066282746310005d0
      b=z+10.5d0
      b=(z+0.5d0)*dlog(b)-b
      do i=1,11
        z=z+1.d0
        a=a+c0(i)/z
      enddo
      lnfacd=b+dlog(cp*a)
      return
      end
 
c  gxurcd0.f to compute Gaunt coefficients a(-m,n,m,v,p)
cu uses lnfacd.f to compute ln(z!)
	subroutine gxurcd0(m,n,v,qmax,na)
	implicit double precision (a-h,o-z)
	include 'gmm01f.par'
	parameter (ng0=np*(2*np**3+10*np**2+19*np+5)/6)
	integer v,qmax,p
	double precision lnfacd,cnv(np,np),ga0(ng0)
	common/cfnv/cnv
	common/g0/ga0
	fb(n,v,p)=dble(p-(n+v+1))*dble(p+(n+v+1))
     +           *dble(p-(n-v))*dble(p+(n-v))
     +           /(dble(2*p+1)*dble(2*p-1))	
	if(iabs(m).gt.n.or.iabs(m).gt.v) then
	   write(6,*) 'warning: |m|>n or v in gxurcd0'
	   qmax=-1	   
	   return
	endif
	qmax=min(n,v)
	nq=qmax+1
	if(n.le.v) then 
	   c1=cnv(n,v)
	else
	   c1=cnv(v,n)
	endif
	c1=c1-lnfacd(dble(n-m))-lnfacd(dble(v+m))
	ga0(na+1)=dexp(c1)
	if(qmax.lt.1) return	
 	p=n+v
	do 8 i=2,nq
	   p=p-2
	   if(m.eq.0) then
	      c1=fb(n,v,p+1)
	      c2=fb(n,v,p+2)
	      goto 2
	   endif
	   c1=fb(n,v,p+1)
	   c2=dble(4*m*m)+fb(n,v,p+2)+fb(n,v,p+3)
	   if(i.eq.2) goto 2
	   c3=-fb(n,v,p+4)
	   goto 4
  2	   ga0(na+i)=c2*ga0(na+i-1)/c1
	   goto 8
  4	   ga0(na+i)=(c2*ga0(na+i-1)+c3*ga0(na+i-2))/c1
  8     continue
	return
	end

	subroutine gid0(nmax,m,n,iv,id)
	nt=nmax*(nmax+1)*(2*nmax+1)/3+nmax*nmax
	ns=max(1,iabs(m))
 	nc0=nmax-iabs(m)
 	id=nc0*(nc0+1)*(2*nc0+1)/6
	if(m) 10,11,12
 10	id=id+(n-ns)*(nc0+1)+iv-ns+1  
	return
 11	id=id+(n-ns)*nmax+iv
	return
 12	id=id+(nmax-n)*(nc0+1)+nmax-iv
	id=nt-id
	return
	end

      function plgndrd(l,mr,x)
      integer l,mr,m,index,i,ll
      double precision x,fact,pll,pmm,pmmp1,somx2
      double precision plgndrd,lnfacd
      m=mr
      index=0
      if(m.lt.0) then
         index=1
         m=-m
      endif
      if(dabs(x).gt.1.d0) then
	 write(6,*) 'n,m,x: ', l,(-2*index+1)*m,x
         pause 'bad arguments in plgndrd'
      end if
      if(m.gt.l) then
	 plgndrd=0.d0
         return
      end if
      pmm=1.d0
      if(m.gt.0) then
        somx2=dsqrt((1.d0-x)*(1.d0+x))
        fact=1.d0
        do 11 i=1,m
          pmm=-pmm*fact*somx2
          fact=fact+2.d0
11      continue
      endif
      if(l.eq.m) then
        plgndrd=pmm
      else
        pmmp1=x*(2*m+1)*pmm
        if(l.eq.m+1) then
          plgndrd=pmmp1
        else
          do 12 ll=m+2,l
            pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m)
            pmm=pmmp1
            pmmp1=pll
12        continue
          plgndrd=pll
        endif
      endif
      plgndrd=-plgndrd
      if(m/2*2.eq.m) plgndrd=-plgndrd
      if(index.gt.0) then
         fact=lnfacd(dble(l-m))-lnfacd(dble(l+m))
         fact=dexp(fact)
         plgndrd=-plgndrd*fact
         if(m/2*2.eq.m) plgndrd=-plgndrd
      endif 
      return
      end
C

      subroutine rtrT(anpt,nodrj,nodri,ekt,drot,ij1,ij2,ii1,ii2)
      implicit double precision (a-h,o-z)
      include 'gmm01f.par'
      parameter(nmp=np*(np+2))
      parameter (nrc=4*np*(np+1)*(np+2)/3+np,nij=nLp*(nLp-1)/2)
      double precision drot(nrc)
      complex*16 anpt(2,nmp),ant(2,2*np),amt(2,-np:np),a,b,
     +           ekt(np),ek(-np:np),atr(2,np,nmp)
      common/tran/atr
      nj1=dsqrt(dble(ij1))
      nj2=dsqrt(dble(ij2))
      ni1=dsqrt(dble(ii1))
      ni2=dsqrt(dble(ii2))
      ek(0)=1.d0
      nmax=max(nj2,ni2)
      do m=1,nmax
         ek(m)=ekt(m)
         ek(-m)=dconjg(ek(m))
      enddo
      irc=0
      do n=1,nj1-1
         do k=-n,n
            do m=-n,n
               irc=irc+1
            enddo
         enddo
      enddo
      do n=nj1,nj2
         n1=n*(n+1)
         do m=-n,n
            amt(1,m)=0.d0
            amt(2,m)=0.d0
         enddo
         do 10 k=-n,n
            kn=n1+k
            if(kn.lt.ij1) then
              irc=irc+2*n+1
              goto 10
            endif
            if(kn.gt.ij2) goto 11
            a=ek(k)*anpt(1,kn)
            b=ek(k)*anpt(2,kn)
            do m=-n,n
               irc=irc+1
               amt(1,m)=amt(1,m)+a*drot(irc)
               amt(2,m)=amt(2,m)+b*drot(irc)
            enddo
 10      continue
 11      do m=-n,n
            imn=n1+m
            anpt(1,imn)=amt(1,m)
            anpt(2,imn)=amt(2,m)
         enddo
      enddo
      mmax=min(nj2,ni2)
      do m=-mmax,mmax
         n1=max(1,iabs(m))
         n1j=max(n1,nj1)
         do n=n1j,nj2
            imn=n*(n+1)+m
            do ip=1,2
               ant(ip,n)=anpt(ip,imn)
            enddo
         enddo
         n1i=max(n1,ni1)
         do n=n1i,ni2
            imn=n*(n+1)+m
            a=0.d0
            b=0.d0
            do 20 l=n1j,nj2
               ml=l*(l+1)+m
               a=a+atr(1,n,ml)*ant(1,l)
     1            +atr(2,n,ml)*ant(2,l)
               b=b+atr(1,n,ml)*ant(2,l)
     1            +atr(2,n,ml)*ant(1,l)
 20         continue
            anpt(1,imn) = a
            anpt(2,imn) = b
         enddo
      enddo
      in=1
      irc=0
      do n=1,ni1-1
         do k=-n,n
            do m=-n,n
               irc=irc+1
            enddo
         enddo
      enddo
      do n=ni1,ni2
         in=-in
         n1=n*(n+1)
         do m=-n,n
            amt(1,m)=0.d0
            amt(2,m)=0.d0
         enddo
         sik=-in
         do 30 k=-n,n
	    sik=-sik
	    kn=n1+k
	    a=sik*anpt(1,kn)
            b=sik*anpt(2,kn)
            do m=-n,n
               irc=irc+1
               amt(1,m)=amt(1,m)+a*drot(irc)
               amt(2,m)=amt(2,m)+b*drot(irc)
            enddo
 30      continue
         sik=-in
         do 31 m=-n,n
            sik=-sik
            imn=n1+m
            if(imn.lt.ii1.or.imn.gt.ii2) goto 31
            anpt(1,imn)=amt(1,m)*ek(-m)*sik
            anpt(2,imn)=amt(2,m)*ek(-m)*sik
 31      continue
      enddo
      return
      end

	subroutine transT(nL,r0,nmax,uvmax,fint,atr0,btr0,ek,
     +             drot,as,bs,as1,bs1,ind,confg,iuvc,isw)
	implicit double precision (a-h,o-z)
	include 'gmm01f.par'
	parameter (nmp=np*(np+2))
	parameter (ni0=np*(np+1)*(2*np+1)/3+np*np)	
	parameter (nrc=4*np*(np+1)*(np+2)/3+np)
	parameter (nij=nLp*(nLp-1)/2)
	integer v,nmax(nLp),uvmax(nLp),ind(nLp)
	double precision r0(6,nLp),drot(nrc,nij),confg(5,nij)
	complex*16 atr(2,np,nmp),atr0(ni0,nij),btr0(ni0,nij),
     +     ek(np,nij),as(nLp,nmp),bs(nLp,nmp),as1(nLp,nmp),
     +     bs1(nLp,nmp),at1(2,nmp)
        common/tran/atr
        do i=1,nL
           do imn=1,uvmax(i)
              as1(i,imn)=dcmplx(0.d0,0.d0)
	      bs1(i,imn)=dcmplx(0.d0,0.d0)
	   enddo
	enddo
 	do 11 i=1,nL-1
 	   if(ind(i).gt.0) goto 11
	   do 10 j=i+1,nL
	      if(ind(j).gt.0) goto 10
	      ij=(j-1)*(j-2)/2+j-i
	      x0=confg(1,ij)
              y0=confg(2,ij)
              z0=confg(3,ij)
              temp=confg(5,ij)
	      if(temp.le.fint) goto 10
              nlarge=max(nmax(i),nmax(j))
              itrc=0
              nsmall=min(nmax(i),nmax(j))
              do m=-nsmall,nsmall
                 n1=max(1,iabs(m))
                 do n=n1,nlarge
                    do v=n1,nlarge
                       itrc=itrc+1
                       iuv=v*(v+1)+m
                       atr(1,n,iuv)=atr0(itrc,ij)
                       atr(2,n,iuv)=btr0(itrc,ij)
                       if(x0.eq.0.d0.and.y0.eq.0.d0) then
                          if(z0.lt.0.d0) goto 20
                       endif
                       goto 21
 20                    sic=dble((-1)**(n+v))
                       atr(1,n,iuv)=sic*atr(1,n,iuv)
                       atr(2,n,iuv)=-sic*atr(2,n,iuv)
 21                    continue
                    enddo
                 enddo
              enddo
              do iuv=1,uvmax(j)
                 at1(1,iuv)=as(j,iuv)
                 at1(2,iuv)=bs(j,iuv)
              enddo
              If(nmax(j).lt.nlarge) then
                 do n=nmax(j)+1,nlarge
                    do m=-n,n
                       iuv=n*(n+1)+m
                       at1(1,iuv)=dcmplx(0.d0,0.d0)
                       at1(2,iuv)=dcmplx(0.d0,0.d0)
                    enddo
                 enddo
              endif
              goto(22,23),isw
 22	      ij1=iuvc
              ij2=uvmax(j)
              ii1=iuvc
              ii2=uvmax(i)
              goto 24
 23	      ij1=1
              ij2=iuvc-1
              ii1=iuvc
              ii2=uvmax(i)
 24           if(x0.eq.0.d0.and.y0.eq.0.d0) then
                 call trvT(at1,nmax(j),nmax(i),ij1,ij2,ii1,ii2)
              else
                 call rtrT(at1,nmax(j),nmax(i),ek(1,ij),
     +                     drot(1,ij),ij1,ij2,ii1,ii2)
              endif
     	      do imn=ii1,ii2
     	         as1(i,imn)=as1(i,imn)+at1(1,imn)
     	         bs1(i,imn)=bs1(i,imn)+at1(2,imn)
     	      enddo
     	      do m=-nsmall,nsmall
                 n1=max(1,iabs(m))
                 do n=n1,nlarge
                    do v=n1,nlarge
                       iuv=v*(v+1)+m
                       sic=dble((-1)**(n+v))
                       atr(1,n,iuv)=sic*atr(1,n,iuv)
                       atr(2,n,iuv)=-sic*atr(2,n,iuv)
                    enddo
                 enddo
              enddo
              do iuv=1,uvmax(i)
                 at1(1,iuv)=as(i,iuv)
                 at1(2,iuv)=bs(i,iuv)
              enddo
              If(nmax(i).lt.nlarge) then
                 do n=nmax(i)+1,nlarge
                    do m=-n,n
                       iuv=n*(n+1)+m
                       at1(1,iuv)=dcmplx(0.d0,0.d0)
                       at1(2,iuv)=dcmplx(0.d0,0.d0)
                    enddo
                 enddo
              endif
              goto(25,26),isw
 25	      ii1=iuvc
              ii2=uvmax(i)
              ij1=iuvc
              ij2=uvmax(j)
              goto 27
 26	      ii1=1
              ii2=iuvc-1
              ij1=iuvc
              ij2=uvmax(j)
 27           if(x0.eq.0.d0.and.y0.eq.0.d0) then
                 call trvT(at1,nmax(i),nmax(j),ii1,ii2,ij1,ij2)
              else 
                 call rtrT(at1,nmax(i),nmax(j),ek(1,ij),
     +                     drot(1,ij),ii1,ii2,ij1,ij2)
              endif
     	      do imn=ij1,ij2
     	         as1(j,imn)=as1(j,imn)+at1(1,imn)
     	         bs1(j,imn)=bs1(j,imn)+at1(2,imn)
     	      enddo
 10	   continue
 11	continue
	return
	end

      subroutine trvT(anpt,nodrj,nodri,ij1,ij2,ii1,ii2)
      implicit double precision (a-h,o-z)
      include 'gmm01f.par'
      parameter(nmp=np*(np+2))
      complex*16 anpt(2,nmp),ant(2,2*np),a,b,atr(2,np,nmp)
      common/tran/atr

      nj1=dsqrt(dble(ij1))
      nj2=dsqrt(dble(ij2))
      ni1=dsqrt(dble(ii1))
      ni2=dsqrt(dble(ii2))
      mmax=min(nj2,ni2)
      do m=-mmax,mmax
         n1=max(1,iabs(m))
         n1j=max(n1,nj1)
         do 10 n=n1j,nj2
            imn=n*(n+1)+m
            if(imn.lt.ij1.or.imn.gt.ij2) goto 10
            do ip=1,2
               ant(ip,n)=anpt(ip,imn)
            enddo
 10      continue
         n1i=max(n1,ni1)
         do 20 n=n1i,ni2
            imn=n*(n+1)+m
            if(imn.lt.ii1.or.imn.gt.ii2) goto 20
            a=0.d0
            b=0.d0
            do 21 l=n1j,nodrj
               ml=l*(l+1)+m
               if(ml.lt.ij1.or.ml.gt.ij2) goto 21               
               a=a+atr(1,n,ml)*ant(1,l)
     +            +atr(2,n,ml)*ant(2,l)
               b=b+atr(1,n,ml)*ant(2,l)
     +            +atr(2,n,ml)*ant(1,l)
 21         continue
            anpt(1,imn) = a
            anpt(2,imn) = b
 20      continue
      enddo
      return
      end


c
c  subroutine xuwigd.f
c  returns the entire set of Wigner 3jm symbols in double 
c  precision for a given integer group (j1,j2,m1,m2) 
c  needs lnfacd.f to compute ln(z!)
c  using the algorithm described in [Xu, Journal of Computational
c  Physics 139, 137-165 (1998)]
c 
      subroutine xuwigd(j1,j2,m1,m2,c,cf,n,kmax)
      include 'gmm01f.par'
      implicit double precision (a-h,o-z)	
      double precision c(n),cf(n),lnfacd
      common/fnr/fnr(0:4*(np+1))
      data small/1.d-12/
      a(j1,j2,j3,m1,m2)=fnr(j3+j1-j2)*fnr(j1+j2+1+j3)
     +                  *fnr(j3+m1+m2)*fnr(j3-j1+j2)
     +                  *fnr(j1+j2+1-j3)*fnr(j3-m1-m2)
      b(j1,j2,j3,m1,m2)=dble(2*j3+1)*(dble((m1+m2)
     +                  *(j1*j1+j1-j2*j2-j2))
     +                  -dble((m1-m2)*(j3*j3+j3)))
      do i=1,n
         c(i)=0.d0
      enddo
      if(iabs(m1).gt.j1.or.iabs(m2).gt.j2) return
      kmax=j1+j2-max(iabs(j1-j2),iabs(m1+m2))+1
      if(kmax.gt.n) stop
      j3=j1+j2
      m3=-m1-m2
      s=lnfacd(dble(2*j1))+lnfacd(dble(2*j2))
      s=s-lnfacd(dble(2*j1+2*j2+1))+lnfacd(dble(j1+j2-m1-m2))
      s=s-lnfacd(dble(j1-m1))-lnfacd(dble(j1+m1))
      s=s+lnfacd(dble(j1+j2+m1+m2))-lnfacd(dble(j2-m2))
      s=s-lnfacd(dble(j2+m2))
      j=j1-j2+m1+m2
      c(1)=((-1)**j)*dexp(.5d0*s)
      if(kmax.eq.1) return
      c(2)=-b(j1,j2,j3,m1,m2)*c(1)
     +      /(dble(j3+1)*a(j1,j2,j3,m1,m2))
      if(kmax.eq.2) return
      j3f=max(iabs(j1-j2),iabs(m1+m2))
      j3=j3f
      if(j3.eq.0) then
         c(kmax)=dble((-1)**(j1-m1))/fnr(2*j1+1)
      else
	 s=lnfacd(dble(j3+j1-j2))+lnfacd(dble(j3-j1+j2))
	 s=s+lnfacd(dble(j1+j2-j3))+lnfacd(dble(j1-m1))
	 s=s+lnfacd(dble(j1+m1))+lnfacd(dble(j2-m2))
	 s=s+lnfacd(dble(j2+m2))-lnfacd(dble(j1+j2+j3+1))
	 s=s+lnfacd(dble(j3-m3))+lnfacd(dble(j3+m3))
	 t=0.5d0*s
	 if(j3.eq.j1-j2.or.j3.eq.m3) k=j2+m2
	 if(j3.eq.j2-j1.or.j3.eq.-m3) k=j1-m1
	 s=lnfacd(dble(j1+j2-j3-k))+lnfacd(dble(j1-m1-k))
	 s=s+lnfacd(dble(k))+lnfacd(dble(j2+m2-k))
	 s=s+lnfacd(dble(j3-j2+m1+k))+lnfacd(dble(j3-j1-m2+k))
	 j=j1-j2-m3
	 c(kmax)=((-1)**(k+j))*dexp(t-s)           
      endif
      if(kmax.eq.3) return
      if(j3.eq.0) then
         c(kmax-1)=c(kmax)*dble(m1)/(fnr(j1)*fnr(j1+1))
      else
         c(kmax-1)=-b(j1,j2,j3,m1,m2)*c(kmax)/(dble(j3)*
     +            a(j1,j2,j3+1,m1,m2))
      endif
      if(kmax.eq.4) return
      cf(kmax)=c(kmax)
      cf(kmax-1)=c(kmax-1)
      j3=j3f
      do i=kmax-2,3,-1
         j3=j3+1
         cf(i)=dble(j3+1)*a(j1,j2,j3,m1,m2)*cf(i+2)
     +                   +b(j1,j2,j3,m1,m2)*cf(i+1)
         cf(i)=-cf(i)/(dble(j3)*a(j1,j2,j3+1,m1,m2))
      enddo
      j3=j1+j2
      do i=3,kmax-2
         j3=j3-1
         c(i)=dble(j3)*a(j1,j2,j3+1,m1,m2)*c(i-2)+
     +        b(j1,j2,j3,m1,m2)*c(i-1)
         c(i)=-c(i)/(dble(j3+1)*a(j1,j2,j3,m1,m2))
         cr=dabs(c(i)-cf(i))
         if(dabs(c(i)).gt.1.d0) cr=cr/dabs(c(i))
         if(cr.lt.small) then
            do j=i,kmax-2
               c(j)=cf(j)
            enddo
            return
         endif
      enddo
      return
      end
