tmd.lp.f

Go to the documentation of this file.

C   This is NOT the original code by the authors listed below!
C   Modified 2009-2013 by Jussi Leinonen (jsleinonen@gmail.com)
C   to be compatible with the Python extension interface.
C   The requirement of non-profit use only has been dropped from
C   this code with the permission of M. I. Mishchenko; thus the
C   license has been made open source compatible.
 
C   New release including the LAPACK matrix inversion procedure.
C   We thank Cory Davis (University of Edinburgh) for pointing
C   out the possibility of replacing the proprietary NAG matrix 
C   inversion routine by the public-domain LAPACK equivalent.
 
C   CALCULATION OF LIGHT SCATTERING BY POLYDISPERSE, RANDOMLY          
C   ORIENTED PARTICLES OF IDENTICAL AXIALLY SYMMETRIC SHAPE      
 
C   This version of the code uses DOUBLE PRECISION variables
C   and must be used along with the accompanying files tmd.par.f      
C   and lpd.f.    
                                                                       
C   Last update 08/06/2005 
                                                                       
C   The code has been developed by Michael Mishchenko at the NASA      
C   Goddard Institute for Space Studies, New York. This research
C   was funded by the NASA Radiation Sciences Program.
                                                                       
C   The code can be used without limitations in any not-for-       
C   profit scientific research.  We only request that in any      
C   publication using the code the source of the code be acknowledged  
C   and relevant references (see below) be made.                                   
                                                                       
C   This version of the code is applicable to spheroids,               
C   Chebyshev particles, and finite circular cylinders.                
                                                                       
C   The computational method is based on the Watermsn's T-matrix      
C   approach and is described in detail in the following papers:           
C                                                                      
C   1.  M. I. Mishchenko, Light scattering by randomly oriented        
C       axially symmetric particles, J. Opt. Soc. Am. A,               
C       vol. 8, 871-882 (1991).                                        
C                                                                      
C   2.  M. I. Mishchenko, Light scattering by size-shape               
C       distributions of randomly oriented axially symmetric           
C       particles of a size comparable to a wavelength,                
C       Appl. Opt., vol. 32, 4652-4666 (1993).                         
C                                                                      
C   3.  M. I. Mishchenko and L. D. Travis, T-matrix computations       
C       of light scattering by large spheroidal particles,             
C       Opt. Commun., vol. 109, 16-21 (1994).                          
C                                                                      
C   4.  M. I. Mishchenko, L. D. Travis, and A. Macke, Scattering       
C       of light by polydisperse, randomly oriented, finite            
C       circular cylinders, Appl. Opt., vol. 35, 4927-4940 (1996).     
C                                                                      
C   5.  D. J. Wielaard, M. I. Mishchenko, A. Macke, and B. E. Carlson, 
C       Improved T-matrix computations for large, nonabsorbing and    
C       weakly absorbing nonspherical particles and comparison         
C       with geometrical optics approximation, Appl. Opt., vol. 36,    
C       4305-4313 (1997).                                                  
C                                                                      
C   A general review of the T-matrix approach can be found in          
C                                                                      
C   6.  M. I. Mishchenko, L. D. Travis, and D. W. Mackowski,           
C       T-matrix computations of light scattering by nonspherical      
C       particles: a review, J. Quant. Spectrosc. Radiat.             
C       Transfer, vol. 55, 535-575 (1996).                             
C                                                                      
C   The following paper provides a detailed user guide to the          
C   T-matrix code:                                                     
C                                                                      
C   7.  M. I. Mishchenko and L. D. Travis, Capabilities and            
C       limitations of a current FORTRAN implementation of the         
C       T-matrix method for randomly oriented, rotationally            
C       symmetric scatterers, J. Quant. Spectrosc. Radiat. Transfer,   
C       vol. 60, 309-324 (1998).                                       
C
C   These papers are available in the .pdf format at the web site 
C
C   http://www.giss.nasa.gov/~crmim/publications/
C
C   or in hardcopy upon request from Michael Mishchenko    
C   Please e-mail your request to crmim@giss.nasa.gov.      
C
C   A comprehensive book "Scattering, Absorption, and Emission of  
C   Light by Small Particles" (Cambridge University Press, Cambridge,
C   2002) is also available in the .pdf format at the web site 
C
C   http://www.giss.nasa.gov/~crmim/books.html
 
C   Analytical averaging over particle orientations (Ref. 1) makes     
C   this method the fastest exact technique currently available.       
C   The use of an automatic convergence procedure                      
C   (Ref. 2) makes the code convenient in massive computations.        
C   Ref. 4 describes features specific for finite cylinders as         
C   particles with sharp rectangular edges.  Ref. 5 describes further      
C   numerical improvements.                                               
                                                                       
C   The use of extended precision variables (Ref. 3) can                     
C   significantly increase the maximal convergent equivalent-sphere
C   size parameter and make it greater than 200 (depending on
C   refractive index and aspect ratio). The extended-precision code
C   is also available. However, the use of extended precision varibales     
C   results in a greater consumption of CPU time.                  
C   On IBM RISC workstations, that code is approximately               
C   five times slower than this double-precision code.  The            
C   CPU time difference between the double-precision and extended-     
C   precision codes can be larger on supercomputers.                   
                                                                       
C   This is the first part of the full T-matrix code.  The second part,    
C   lpd.f, is completely independent of the firsti part. It contains no    
C   T-matrix-specific subroutines and can be compiled separately.
C   The second part of the code replaces the previously implemented 
C   standard matrix inversion scheme based on Gaussian elimination 
C   by a scheme based on the LU factorization technique.     
C   As described in Ref. 5 above, the use of the LU factorization is          
C   especially beneficial for nonabsorbing or weakly absorbing particles.     
C   In this code we use the LAPACK implementation of the LU factorization
C   scheme. LAPACK stands for Linear Algebra PACKage. The latter is 
C   publicly available at the following internet site:
C
C   http://www.netlib.org/lapack/      
                                                                      
C   INPUT PARAMETERS:                                                  
C                                                                      
C      RAT = 1 - particle size is specified in terms of the            
C                equal-volume-sphere radius                             
C      RAT.NE.1 - particle size is specified in terms of the           
C                equal-surface-area-sphere radius                      
C      NDISTR specifies the distribution of equivalent-sphere radii   
C      NDISTR = 1 - modified gamma distribution                        
C           [Eq. (40) of Ref. 7]                                       
C               AXI=alpha                                              
C               B=r_c                                                  
C               GAM=gamma                                              
C      NDISTR = 2 - log-normal distribution                            
C           [Eq. 41) of Ref. 7]                                        
C               AXI=r_g                                                
C               B=[ln(sigma_g)]**2                                    
C      NDISTR = 3 - power law distribution                             
C           [Eq. (42) of Ref. 7]                                       
C                AXI=r_eff (effective radius)                   
C                B=v_eff (effective variance)                
C                Parameters R1 and R2 (see below) are calculated       
C                automatically for given AXI and B
C      NDISTR = 4 - gamma distribution                                 
C           [Eq. (39) of Ref. 7]                                       
C                AXI=a                                                 
C                B=b                                                   
C      NDISTR = 5 - modified power law distribution
C         [Eq. (24) in M. I. Mishchenko et al.,
C         Bidirectional reflectance of flat,
C         optically thick particulate laters: an efficient radiative
C         transfer solution and applications to snow and soil surfaces,
C         J. Quant. Spectrosc. Radiat. Transfer, Vol. 63, 409-432 (1999)].
C                B=alpha
C                                                                      
C      The code computes NPNAX size distributions of the same type     
C      and with the same values of B and GAM in one run.               
C      The parameter AXI varies from AXMAX to AXMAX/NPNAX in steps of  
C      AXMAX/NPNAX.  To compute a single size distribution, use        
C      NPNAX=1 and AXMAX equal to AXI of this size distribution.       
C                                                                      
C      R1 and R2 - minimum and maximum equivalent-sphere radii         
C           in the size distribution. They are calculated automatically              
C           for the power law distribution with given AXI and B   
C           but must be specified for other distributions        
C           after the lines           
C                                                                      
C             DO 600 IAX=1,NPNAX                            
C                AXI=AXMAX-DAX*DFLOAT(IAX-1)                
C                                                                      
C           in the main program.                                 
C           For the modified power law distribution (NDISTR=5), the
C           minimum radius is 0, R2 is the maximum radius,
C           and R1 is the intermediate radius at which the
C           n(r)=const dependence is replaced by the power law
C           dependence.
C                                                                      
C      NKMAX.LE.988 is such that NKMAX+2 is the                        
C           number of Gaussian quadrature points used in               
C           integrating over the size distribution for particles with
C           AXI=AXMAX.  For particles with AXI=AXMAX-AXMAX/NPNAX,      
C           AXMAX-2*AXMAX/NPNAX, etc. the number of Gaussian points    
C           linearly decreases.                                       
C           For the modified power law distribution, the number
C           of integration points on the interval [0,R1] is also
C           equal to NKMAX.
C                                                                      
C      LAM - wavelength of light                                       
C      MRR and MRI - real and imaginary parts of the refractive        
C                  index (MRI.GE.0)   
C      EPS and NP - specify the shape of the particles.                
C             For spheroids NP=-1 and EPS is the ratio of the          
C                 horizontal to rotational axes.  EPS is larger than   
C                 1 for oblate spheroids and smaller than 1 for       
C                 prolate spheroids.                                   
C             For cylinders NP=-2 and EPS is the ratio of the          
C                 diameter to the length.                              
C             For Chebyshev particles NP must be positive and 
C                 is the degree of the Chebyshev polynomial, while     
C                 EPS is the deformation parameter                     
C                 [Eq. (33) of Ref. 7].      
C      DDELT - accuracy of the computations                            
C      NPNA - number of equidistant scattering angles (from 0      
C             to 180 deg) for which the scattering matrix is           
C             calculated.                                              
C      NDGS - parameter controlling the number of division points      
C             in computing integrals over the particle surface.        
C             For compact particles, the recommended value is 2.       
C             For highly aspherical particles larger values (3, 4,...) 
C             may be necessary to obtain convergence.                  
C             The code does not check convergence over this parameter. 
C             Therefore, control comparisons of results obtained with  
C             different NDGS-values are recommended.                   
                                                                       
 
C   OUTPUT PARAMETERS:                                                 
C                                                                      
C      REFF and VEFF - effective radius and effective variance of      
C          the size distribution as defined by Eqs. (43)-(45) of
C          Ref. 7.       
C      CEXT - extinction cross section per particle                    
C      CSCA - scattering cross section per particle                    
C      W - single scattering albedo                                    
C      <cos> - asymmetry parameter of the phase function               
C      ALPHA1,...,BETA2 - coefficients appearing in the expansions
C          of the elements of the scattering matrix in
C          generalized spherical functions
C          [Eqs. (11)-(16) of Ref. 7].
C      F11,...,F44 - elements of the normalized scattering matrix [as      
C          defined by Eqs. (1)-(3) of Ref. 7] versus scattering angle
 
C   Note that LAM, r_c, r_g, r_eff, a, R1, and R2 must 
C   be given in the same units of length, and that 
C   the dimension of CEXT and CSCA is that of LAM squared (e.g., if        
C   LAM and AXI are given in microns, then CEXT and CSCA are         
C   calculated in square microns).    
                                                                       
C   The physical correctness of the computed results is tested using   
C   the general inequalities derived by van der Mee and Hovenier,      
C   Astron. Astrophys., vol. 228, 559-568 (1990).  Although            
C   the message that the test of van der Mee and Hovenier is satisfied 
C   does not guarantee that the results are absolutely correct,        
C   the message that the test is not satisfied can mean that something 
C   is wrong.                                                          
                                                                       
C   The convergence of the T-matrix method for particles with          
C   different sizes, refractive indices, and aspect ratios can be      
C   dramatically different.  Usually, large sizes and large aspect     
C   ratios cause problems.  The user of this code                      
C   should first experiment with different input parameters in          
C   order to get an idea of the range of applicability of this         
C   technique.  Sometimes decreasing the aspect ratio                  
C   from 3 to 2 can increase the maximum convergent equivalent-        
C   sphere size parameter by a factor of several (Ref. 7).                      
C   The CPU time required rapidly increases with increasing ratio      
C   radius/wavelength and/or with increasing particle asphericity.     
C   This should be taken into account in planning massive computations.
C   Using an optimizing compiler on IBM RISC workstations saves        
C   about 70% of CPU time.                                             
                                                                       
C   Execution can be automatically terminated if dimensions of certain 
C   arrays are not big enough or if the convergence procedure decides  
C   that the accuracy of double precision variables is insufficient    
C   to obtain a converged T-matrix solution for given particles.       
C   In all cases, a message appears explaining the cause of termination. 
                                                                       
C   The message                                                        
C        "WARNING:  W IS GREATER THAN 1"                               
C   means that the single-scattering albedo exceeds the maximum        
C   possible value 1.  If W is greater than 1 by more than             
C   DDELT, this message can be an indication of numerical              
C   instability caused by extreme values of particle parameters.       
                                                                       
C   The message "WARNING: NGAUSS=NPNG1" means that convergence over    
C   the parameter NG (see Ref. 2) cannot be obtained for the NPNG1     
C   value specified in the PARAMETER statement in the file tmd.par.f. 
C   Often this is not a serious problem, especially for compact         
C   particles.                                                         
                                                                       
C   Larger and/or more aspherical particles may require larger         
C   values of the parameters NPN1, NPN4, and NPNG1 in the file
C   tmd.par.f.  It is recommended to keep NPN1=NPN4+25 and
C   NPNG1=3*NPN1.  Note that the memory requirement increases    
C   as the third power of NPN4. If the memory of                  
C   a computer is too small to accomodate the code in its current    
C   setting, the parameters NPN1, NPN4, and NPNG1 should be
C   decreased. However, this will decrease the maximum size parameter  
C   that can be handled by the code.                                   
                                                                       
C   In some cases any increases of NPN1 will not make the T-matrix     
C   computations convergent.  This means that the particle is just     
C   too "bad" (extreme size parameter and/or extreme aspect ratio      
C   and/or extreme refractive index; see Ref. 7).              
C   The main program contains several PRINT statements which are       
C   currently commentd out.  If uncommented, these statements will     
C   produce numbers which show the convergence rate and can be         
C   used to determine whether T-matrix computations for given particle 
C   parameters will converge at all.   
                                                                       
C   Some of the common blocks are used to save memory rather than      
C   to transfer data.  Therefore, if a compiler produces a warning     
C   message that the lengths of a common block are different in        
C   different subroutines, this is not a real problem.                 
                                                                       
C   The recommended value of DDELT is 0.001.  For bigger values,       
C   false convergence can be obtained.                                 
                                                                       
C   In computations for spheres use EPS=1.000001 instead of EPS=1.     
C   The use of EPS=1 can cause overflows in some rare cases.           
                                                                       
C   To calculate a monodisperse particle, use the options              
C        NPNAX=1                                                       
C        AXMAX=R                                                       
C        B=1D-1   
C        NKMAX=-1                                                      
C        NDISTR=4                                                      
C        ...                                                           
C        DO 600 IAX=1,NPNAX                                            
C           AXI=AXMAX-DAX*DFLOAT(IAX-1)                                
C           R1=0.9999999*AXI                                           
C           R2=1.0000001*AXI                                          
C        ...                                                           
C   where R is the equivalent-sphere radius.                           
                                                                       
C   It is recommended to use the power law rather than the             
C   gamma size distribution, because in this case convergent solution  
C   can be obtained for larger REFF and VEFF assuming the same
C   maximal R2 (Mishchenko and Travis, Appl. Opt., vol. 33, 7206-7225,
C   1994).           
                                                                       
C   For some compilers, DACOS must be raplaced by DARCOS and DASIN     
C   by DARSIN.                                                         
                                                                       
C   If many different size distributions are computed and the          
C   refractive index is fixed, then another approach can be more       
C   efficient than running this code many times.  Specifically,        
C   scattering results should be computed for monodisperse particles   
C   with sizes ranging from essentially zero to some maximum value     
C   with a small step size (say, 0.02 microns).  These results         
C   should be stored on disk and can be used along with spline         
C   interpolation to compute scattering by particles with intermediate 
C   sizes.  Scattering patterns for monodisperse nonspherical          
C   particles in random orientation are (much) smoother than for 
C   monodisperse spheres, and spline interpolation usually gives good 
C   results. In this way, averaging over any size distribution is a
C   matter of seconds.  For more on size averaging, see Refs. 2 and 4. 
                                                                       
C   We would highly appreciate informing me of any problems encountered 
C   with this code.  Please send your message to the following         
C   e-mail address:  CRMIM@GISS.NASA.GOV.                              
 
C   WHILE THE COMPUTER PROGRAM HAS BEEN TESTED FOR A VARIETY OF CASES,
C   IT IS NOT INCONCEIVABLE THAT IT CONTAINS UNDETECTED ERRORS. ALSO,
C   INPUT PARAMETERS CAN BE USED WHICH ARE OUTSIDE THE ENVELOPE OF
C   VALUES FOR WHICH RESULTS ARE COMPUTED ACCURATELY. FOR THIS REASON,
C   THE AUTHORS AND THEIR ORGANIZATION DISCLAIM ALL LIABILITY FOR
C   ANY DAMAGES THAT MAY RESULT FROM THE USE OF THE PROGRAM. 
                                                                       
      SUBROUTINE calcrand(RAT,LAM,MRR,MRI,EPS,NP,DDELT,NDGS,
     &     NPNAX,AXMAX,B,GAM,NDISTR,NKMAX,NPNA,NCOEFF,
     &     REFF,VEFF,CEXTIN,CSCAT,WALB,ASYMM, F)
      IMPLICIT REAL*8 (a-h,o-z)
      include 'tmd.par.f'
      real*8 reff,veff,cextin,cscat,walb,asymm
      real*8 alpha(npl,4),beta(npl,2),f(npna,4,4)
Cf2py intent(in) rat
Cf2py intent(in) lam
Cf2py intent(in) mrr
Cf2py intent(in) mri
Cf2py intent(in) eps
Cf2py intent(in) np
Cf2py intent(in) ddelt
Cf2py intent(in) ndgs
Cf2py intent(in) npnax
Cf2py intent(in) axmax
Cf2py intent(in) b
Cf2py intent(in) gam
Cf2py intent(in) ndistr
Cf2py intent(in) nkmax
Cf2py intent(in) npna
Cf2py intent(in) ncoeff
Cf2py intent(out) reff
Cf2py intent(out) veff
Cf2py intent(out) cextin
Cf2py intent(out) cscat
Cf2py intent(out) walb
Cf2py intent(out) asymm
Cf2py intent(out) alpha
Cf2py intent(out) beta
Cf2py intent(out) f
 
      real*8  lam,mrr,mri,x(npng2),w(npng2),s(npng2),ss(npng2),
     *        an(npn1),r(npng2),dr(npng2),
     *        ddr(npng2),drr(npng2),dri(npng2),ann(npn1,npn1)
      real*8  xg(1000),wg(1000),tr1(npn2,npn2),ti1(npn2,npn2),
     &        alph1(npl),alph2(npl),alph3(npl),alph4(npl),bet1(npl),
     &        bet2(npl),xg1(2000),wg1(2000),
     &        al1(npl),al2(npl),al3(npl),al4(npl),be1(npl),be2(npl)
      real*4
     &     rt11(npn6,npn4,npn4),rt12(npn6,npn4,npn4),
     &     rt21(npn6,npn4,npn4),rt22(npn6,npn4,npn4),
     &     it11(npn6,npn4,npn4),it12(npn6,npn4,npn4),
     &     it21(npn6,npn4,npn4),it22(npn6,npn4,npn4)
 
 
      COMMON /ct/ tr1,ti1
      COMMON /tmat/ rt11,rt12,rt21,rt22,it11,it12,it21,it22
      p=dacos(-1d0)
 
C  OPEN FILES *******************************************************
 
      OPEN (6,file='tmatr.test')
      OPEN (10,file='tmatr.write')
 
C  INPUT DATA ********************************************************
 
C      RAT=0.5 D0
C      NDISTR=3
C      AXMAX=1D0
C      NPNAX=2
C      B=0.1D0
C      GAM=0.5D0
C      NKMAX=5
C      EPS=2D0
C      NP=-1
C      LAM=0.5D0
C      MRR=1.53 d0
C      MRI=0.008D0
C      DDELT=0.001D0
C      NPNA=19
C      NDGS=2
 
      ncheck=0
      IF (np.EQ.-1.OR.np.EQ.-2) ncheck=1
      IF (np.GT.0.AND.(-1)**np.EQ.1) ncheck=1
      WRITE (6,5454) ncheck
 5454 FORMAT ('NCHECK=',i1)
      dax=axmax/npnax
      IF (abs(rat-1d0).GT.1d-8.AND.np.EQ.-1) CALL sarea (eps,rat)
      if (abs(rat-1d0).GT.1d-8.AND.np.GE.0) CALL surfch(np,eps,rat)
      IF (abs(rat-1d0).GT.1d-8.AND.np.EQ.-2) CALL sareac (eps,rat)
C     PRINT 8000, RAT
 8000 FORMAT ('RAT=',f8.6)
      IF(np.EQ.-1.AND.eps.GE.1d0) print 7000,eps
      IF(np.EQ.-1.AND.eps.LT.1d0) print 7001,eps
      IF(np.GE.0) print 7100,np,eps
      IF(np.EQ.-2.AND.eps.GE.1d0) print 7150,eps
      IF(np.EQ.-2.AND.eps.LT.1d0) print 7151,eps
      print 7400, lam,mrr,mri
      print 7200, ddelt
 7000 FORMAT('RANDOMLY ORIENTED OBLATE SPHEROIDS, A/B=',f11.7)
 7001 FORMAT('RANDOMLY ORIENTED PROLATE SPHEROIDS, A/B=',f11.7)
 7100 FORMAT('RANDOMLY ORIENTED CHEBYSHEV PARTICLES, T',
     &       i1,'(',f5.2,')')
 7150 FORMAT('RANDOMLY ORIENTED OBLATE CYLINDERS, D/L=',f11.7)
 7151 FORMAT('RANDOMLY ORIENTED PROLATE CYLINDERS, D/L=',f11.7)
 7200 FORMAT ('ACCURACY OF COMPUTATIONS DDELT = ',d8.2)
 7400 FORMAT('LAM=',f10.6,3x,'MRR=',d10.4,3x,'MRI=',d10.4)
      ddelt=0.1d0*ddelt
      DO 600 iax=1,npnax
         axi=axmax-dax*dfloat(iax-1)
         r1=0.89031d0*axi
         r2=1.56538d0*axi
         nk=int(axi*nkmax/axmax+2)
         IF (nk.GT.1000) print 8001,nk
         IF (nk.GT.1000) stop
         IF (ndistr.EQ.3) CALL power (axi,b,r1,r2)
 8001    FORMAT ('NK=',i4,' I.E., IS GREATER THAN 1000. ',
     &           'EXECUTION TERMINATED.')
         CALL gauss (nk,0,0,xg,wg)
         z1=(r2-r1)*0.5d0
         z2=(r1+r2)*0.5d0
         z3=r1*0.5d0
         IF (ndistr.EQ.5) GO TO 3
         DO i=1,nk
            xg1(i)=z1*xg(i)+z2
            wg1(i)=wg(i)*z1
         ENDDO
         GO TO 4
    3    DO i=1,nk
            xg1(i)=z3*xg(i)+z3
            wg1(i)=wg(i)*z3
         ENDDO
         DO i=nk+1,2*nk
            ii=i-nk
            xg1(i)=z1*xg(ii)+z2
            wg1(i)=wg(ii)*z1
         ENDDO
         nk=nk*2
    4    CALL distrb (nk,xg1,wg1,ndistr,axi,b,gam,r1,r2,
     &                reff,veff,p)
         print 8002,r1,r2
 8002    FORMAT('R1=',f10.6,'   R2=',f10.6)
         IF (abs(rat-1d0).LE.1d-6) print 8003, reff,veff
         IF (abs(rat-1d0).GT.1d-6) print 8004, reff,veff
 8003    FORMAT('EQUAL-VOLUME-SPHERE REFF=',f8.4,'   VEFF=',f7.4)
 8004    FORMAT('EQUAL-SURFACE-AREA-SPHERE REFF=',f8.4,
     &          '   VEFF=',f7.4)
         print 7250,nk
 7250    FORMAT('NUMBER OF GAUSSIAN QUADRATURE POINTS ',
     &          'IN SIZE AVERAGING =',i4)
         DO i=1,npl
            alph1(i)=0d0
            alph2(i)=0d0
            alph3(i)=0d0
            alph4(i)=0d0
            bet1(i)=0d0
            bet2(i)=0d0
            alpha(i,1)=0d0
            alpha(i,2)=0d0
            alpha(i,3)=0d0
            alpha(i,4)=0d0
            beta(i,1)=0d0
            beta(i,2)=0d0
         ENDDO      
         cscat=0d0
         cextin=0d0
         l1max=0
         DO 500 ink=1,nk
            i=nk-ink+1
            a=rat*xg1(i)
            xev=2d0*p*a/lam
            ixxx=xev+4.05d0*xev**0.333333d0
            inm1=max0(4,ixxx)
            IF (inm1.GE.npn1) print 7333, npn1
            IF (inm1.GE.npn1) stop
 7333 FORMAT('CONVERGENCE IS NOT OBTAINED FOR NPN1=',i3,  
     &       '.  EXECUTION TERMINATED')
            qext1=0d0
            qsca1=0d0
            DO 50 nma=inm1,npn1
               nmax=nma
               mmax=1
               ngauss=nmax*ndgs
               IF (ngauss.GT.npng1) print 7340, ngauss
               IF (ngauss.GT.npng1) stop
 7340          FORMAT('NGAUSS =',i3,' I.E. IS GREATER THAN NPNG1.',
     &                '  EXECUTION TERMINATED')
 7334          FORMAT(' NMAX =', i3,'  DC2=',d8.2,'   DC1=',d8.2)
 7335 FORMAT('                              NMAX1 =', i3,'  DC2=',d8.2,
     &      '  DC1=',d8.2)
               CALL const(ngauss,nmax,mmax,p,x,w,an,ann,s,ss,np,eps)
               CALL vary(lam,mrr,mri,a,eps,np,ngauss,x,p,ppi,pir,pii,r,
     &                   dr,ddr,drr,dri,nmax)
               CALL tmatr0 (ngauss,x,w,an,ann,s,ss,ppi,pir,pii,r,dr,
     &                      ddr,drr,dri,nmax,ncheck)
               qext=0d0
               qsca=0d0
               DO n=1,nmax
                  n1=n+nmax
                  tr1nn=tr1(n,n)
                  ti1nn=ti1(n,n)
                  tr1nn1=tr1(n1,n1)
                  ti1nn1=ti1(n1,n1)
                  dn1=dfloat(2*n+1)
                  qsca=qsca+dn1*(tr1nn*tr1nn+ti1nn*ti1nn
     &                          +tr1nn1*tr1nn1+ti1nn1*ti1nn1)
                  qext=qext+(tr1nn+tr1nn1)*dn1
               ENDDO
               dsca=abs((qsca1-qsca)/qsca)
               dext=abs((qext1-qext)/qext)
C              PRINT 7334, NMAX,DSCA,DEXT
               qext1=qext
               qsca1=qsca
               nmin=dfloat(nmax)/2d0+1d0
               DO 10 n=nmin,nmax
                  n1=n+nmax
                  tr1nn=tr1(n,n)
                  ti1nn=ti1(n,n)
                  tr1nn1=tr1(n1,n1)
                  ti1nn1=ti1(n1,n1)
                  dn1=dfloat(2*n+1)
                  dqsca=dn1*(tr1nn*tr1nn+ti1nn*ti1nn
     &                      +tr1nn1*tr1nn1+ti1nn1*ti1nn1)
                  dqext=(tr1nn+tr1nn1)*dn1
                  dqsca=abs(dqsca/qsca)
                  dqext=abs(dqext/qext)
                  nmax1=n
                  IF (dqsca.LE.ddelt.AND.dqext.LE.ddelt) GO TO 12
   10          CONTINUE
   12          CONTINUE
c              PRINT 7335, NMAX1,DQSCA,DQEXT
               IF(dsca.LE.ddelt.AND.dext.LE.ddelt) GO TO 55
               IF (nma.EQ.npn1) print 7333, npn1
               IF (nma.EQ.npn1) stop      
   50       CONTINUE
   55       nnnggg=ngauss+1
            IF (ngauss.EQ.npng1) print 7336
            mmax=nmax1
            DO 150 ngaus=nnnggg,npng1
               ngauss=ngaus
               nggg=2*ngauss
 7336          FORMAT('WARNING: NGAUSS=NPNG1')
 7337          FORMAT(' NG=',i3,'  DC2=',d8.2,'   DC1=',d8.2)
               CALL const(ngauss,nmax,mmax,p,x,w,an,ann,s,ss,np,eps)
               CALL vary(lam,mrr,mri,a,eps,np,ngauss,x,p,ppi,pir,pii,r,
     &                   dr,ddr,drr,dri,nmax)
               CALL tmatr0 (ngauss,x,w,an,ann,s,ss,ppi,pir,pii,r,dr,
     &                      ddr,drr,dri,nmax,ncheck)
               qext=0d0
               qsca=0d0
               DO 104 n=1,nmax
                  n1=n+nmax
                  tr1nn=tr1(n,n)
                  ti1nn=ti1(n,n)
                  tr1nn1=tr1(n1,n1)
                  ti1nn1=ti1(n1,n1)
                  dn1=dfloat(2*n+1)
                  qsca=qsca+dn1*(tr1nn*tr1nn+ti1nn*ti1nn
     &                          +tr1nn1*tr1nn1+ti1nn1*ti1nn1)
                  qext=qext+(tr1nn+tr1nn1)*dn1
  104          CONTINUE
               dsca=abs((qsca1-qsca)/qsca)
               dext=abs((qext1-qext)/qext)
c              PRINT 7337, NGGG,DSCA,DEXT
               qext1=qext
               qsca1=qsca
               IF(dsca.LE.ddelt.AND.dext.LE.ddelt) GO TO 155
               IF (ngaus.EQ.npng1) print 7336
  150       CONTINUE
  155       CONTINUE
            qsca=0d0
            qext=0d0
            nnm=nmax*2
            DO 204 n=1,nnm
               qext=qext+tr1(n,n)
  204       CONTINUE
            IF (nmax1.GT.npn4) print 7550, nmax1
 7550       FORMAT ('NMAX1 = ',i3, ', i.e. greater than NPN4.',
     &              ' Execution terminated')
            IF (nmax1.GT.npn4) stop              
            DO 213 n2=1,nmax1
               nn2=n2+nmax
               DO 213 n1=1,nmax1
                  nn1=n1+nmax
                  zz1=tr1(n1,n2)
                  rt11(1,n1,n2)=zz1
                  zz2=ti1(n1,n2)
                  it11(1,n1,n2)=zz2
                  zz3=tr1(n1,nn2)
                  rt12(1,n1,n2)=zz3
                  zz4=ti1(n1,nn2)
                  it12(1,n1,n2)=zz4
                  zz5=tr1(nn1,n2)
                  rt21(1,n1,n2)=zz5
                  zz6=ti1(nn1,n2)
                  it21(1,n1,n2)=zz6
                  zz7=tr1(nn1,nn2)
                  rt22(1,n1,n2)=zz7
                  zz8=ti1(nn1,nn2)
                  it22(1,n1,n2)=zz8
                  qsca=qsca+zz1*zz1+zz2*zz2+zz3*zz3+zz4*zz4
     &                 +zz5*zz5+zz6*zz6+zz7*zz7+zz8*zz8
  213       CONTINUE
C           PRINT 7800,0,ABS(QEXT),QSCA,NMAX
            DO 220 m=1,nmax1
               CALL tmatr(m,ngauss,x,w,an,ann,s,ss,ppi,pir,pii,r,dr,
     &                    ddr,drr,dri,nmax,ncheck)
               nm=nmax-m+1
               nm1=nmax1-m+1
               m1=m+1
               qsc=0d0
               DO 214 n2=1,nm1
                  nn2=n2+m-1
                  n22=n2+nm
                  DO 214 n1=1,nm1
                     nn1=n1+m-1
                     n11=n1+nm
                     zz1=tr1(n1,n2)
                     rt11(m1,nn1,nn2)=zz1
                     zz2=ti1(n1,n2)
                     it11(m1,nn1,nn2)=zz2
                     zz3=tr1(n1,n22)
                     rt12(m1,nn1,nn2)=zz3
                     zz4=ti1(n1,n22)
                     it12(m1,nn1,nn2)=zz4
                     zz5=tr1(n11,n2)
                     rt21(m1,nn1,nn2)=zz5
                     zz6=ti1(n11,n2)
                     it21(m1,nn1,nn2)=zz6
                     zz7=tr1(n11,n22)
                     rt22(m1,nn1,nn2)=zz7
                     zz8=ti1(n11,n22)
                     it22(m1,nn1,nn2)=zz8
                     qsc=qsc+(zz1*zz1+zz2*zz2+zz3*zz3+zz4*zz4
     &                       +zz5*zz5+zz6*zz6+zz7*zz7+zz8*zz8)*2d0
  214          CONTINUE
               nnm=2*nm
               qxt=0d0
               DO 215 n=1,nnm
                  qxt=qxt+tr1(n,n)*2d0
  215          CONTINUE
               qsca=qsca+qsc
               qext=qext+qxt
C              PRINT 7800,M,ABS(QXT),QSC,NMAX
 7800          FORMAT(' m=',i3,'  qxt=',d12.6,'  qsc=',d12.6,
     &                '  nmax=',i3)
  220       CONTINUE
            coeff1=lam*lam*0.5d0/p
            csca=qsca*coeff1
            cext=-qext*coeff1
c           PRINT 7880, NMAX,NMAX1
 7880       FORMAT ('nmax=',i3,'   nmax1=',i3)
            CALL gsp (nmax1,csca,lam,al1,al2,al3,al4,be1,be2,lmax)
            l1m=lmax+1
            l1max=max(l1max,l1m)
            wgii=wg1(i)
            wgi=wgii*csca
            DO 250 l1=1,l1m
               alph1(l1)=alph1(l1)+al1(l1)*wgi
               alph2(l1)=alph2(l1)+al2(l1)*wgi
               alph3(l1)=alph3(l1)+al3(l1)*wgi
               alph4(l1)=alph4(l1)+al4(l1)*wgi
               bet1(l1)=bet1(l1)+be1(l1)*wgi
               bet2(l1)=bet2(l1)+be2(l1)*wgi
  250       CONTINUE
            cscat=cscat+wgi
            cextin=cextin+cext*wgii
C           PRINT 6070, I,NMAX,NMAX1,NGAUSS
 6070       FORMAT(4i6)
  500    CONTINUE
         DO 510 l1=1,l1max
            alph1(l1)=alph1(l1)/cscat
            alph2(l1)=alph2(l1)/cscat
            alph3(l1)=alph3(l1)/cscat
            alph4(l1)=alph4(l1)/cscat
            bet1(l1)=bet1(l1)/cscat
            bet2(l1)=bet2(l1)/cscat
            alpha(l1,1)=alph1(l1)
            alpha(l1,2)=alph2(l1)
            alpha(l1,3)=alph3(l1)
            alpha(l1,4)=alph4(l1)
            beta(l1,1)=bet1(l1)
            beta(l1,2)=bet2(l1)
  510    CONTINUE
         walb=cscat/cextin
         CALL hovenr(l1max,alph1,alph2,alph3,alph4,bet1,bet2)
         asymm=alph1(2)/3d0
         print 9100,cextin,cscat,walb,asymm
 9100    FORMAT('CEXT=',d12.6,2x,'CSCA=',d12.6,2x,
     &          2x,'W=',d12.6,2x,'<COS>=',d12.6)
         IF (walb.GT.1d0) print 9111
 9111    FORMAT ('WARNING: W IS GREATER THAN 1')
         WRITE (10,580) walb,l1max
         DO l=1,l1max
            WRITE (10,575) alph1(l),alph2(l),alph3(l),alph4(l),
     &                     bet1(l),bet2(l)
         ENDDO   
  575    FORMAT(6d14.7)
  580    FORMAT(d14.8,i8)
         lmax=l1max-1
         CALL matr (alph1,alph2,alph3,alph4,bet1,bet2,lmax,npna,f)
  600 CONTINUE
      itime=mclock()
      time=dfloat(itime)/6000d0
      print 1001,time
 1001 FORMAT (' time =',f8.2,' min')
      RETURN
      END
 
C**********************************************************************
 
C**********************************************************************
 
C**********************************************************************
 
C**********************************************************************
 
C**********************************************************************
 
C**********************************************************************
 
C**********************************************************************
 
C**********************************************************************
 
C**********************************************************************
C                                                                     *
C   CALCULATION OF SPHERICAL BESSEL FUNCTIONS OF THE FIRST KIND       *
C   J=JR+I*JI OF COMPLEX ARGUMENT X=XR+I*XI OF ORDERS FROM 1 TO NMAX  *
C   BY USING BACKWARD RECURSION. PARAMETR NNMAX DETERMINES NUMERICAL  *
C   ACCURACY. U=UR+I*UI - FUNCTION (1/X)(D/DX)(X*J(X))                *
C                                                                     *
C**********************************************************************
 
C**********************************************************************
 
C**********************************************************************
 
C**********************************************************************
 
C********************************************************************
C                                                                   *
C   CALCULATION OF THE EXPANSION COEFFICIENTS FOR (I,Q,U,V) -       *
C   REPRESENTATION.                                                 *
C                                                                   *
C   INPUT PARAMETERS:                                               *
C                                                                   *
C      LAM - WAVELENGTH OF LIGHT                                    *
C      CSCA - SCATTERING CROSS SECTION                              *
C      TR AND TI - ELEMENTS OF THE T-MATRIX. TRANSFERRED THROUGH    *
C                  COMMON /CTM/                                     *
C      NMAX - DIMENSION OF T(M)-MATRICES                            *
C                                                                   *
C   OUTPUT INFORTMATION:                                            *
C                                                                   *
C      ALF1,...,ALF4,BET1,BET2 - EXPANSION COEFFICIENTS             *
C      LMAX - NUMBER OF COEFFICIENTS MINUS 1                        *
C                                                                   *
C********************************************************************
 
      SUBROUTINE gsp(NMAX,CSCA,LAM,ALF1,ALF2,ALF3,ALF4,BET1,BET2,LMAX)
      include 'tmd.par.f'
      IMPLICIT REAL*8 (a-b,d-h,o-z),COMPLEX*16 (C)
      real*8 lam,ssign(900)
      real*8  csca,ssi(npl),ssj(npn1),
     &        alf1(npl),alf2(npl),alf3(npl),
     &        alf4(npl),bet1(npl),bet2(npl),
     &        tr1(npl1,npn4),tr2(npl1,npn4),
     &        ti1(npl1,npn4),ti2(npl1,npn4),
     &        g1(npl1,npn6),g2(npl1,npn6),
     &        ar1(npn4),ar2(npn4),ai1(npn4),ai2(npn4),
     &        fr(npn4,npn4),fi(npn4,npn4),ff(npn4,npn4)
      real*4 b1r(npl1,npl1,npn4),b1i(npl1,npl1,npn4),
     &       b2r(npl1,npl1,npn4),b2i(npl1,npl1,npn4),
     &       d1(npl1,npn4,npn4),d2(npl1,npn4,npn4),
     &       d3(npl1,npn4,npn4),d4(npl1,npn4,npn4),
     &       d5r(npl1,npn4,npn4),d5i(npl1,npn4,npn4),
     &       plus1(npn6*npn4*npn4*8)         
      real*4
     &     tr11(npn6,npn4,npn4),tr12(npn6,npn4,npn4),
     &     tr21(npn6,npn4,npn4),tr22(npn6,npn4,npn4),
     &     ti11(npn6,npn4,npn4),ti12(npn6,npn4,npn4),
     &     ti21(npn6,npn4,npn4),ti22(npn6,npn4,npn4)
      COMPLEX*16 CIM(NPN1)
 
      COMMON /tmat/ tr11,tr12,tr21,tr22,ti11,ti12,ti21,ti22
      COMMON /cbess/ b1r,b1i,b2r,b2i    
      COMMON /ss/ ssign
      equivalence( plus1(1),tr11(1,1,1) )
      equivalence(d1(1,1,1),plus1(1)),        
     &            (d2(1,1,1),plus1(npl1*npn4*npn4+1)),
     &            (d3(1,1,1),plus1(npl1*npn4*npn4*2+1)),
     &            (d4(1,1,1),plus1(npl1*npn4*npn4*3+1)), 
     &            (d5r(1,1,1),plus1(npl1*npn4*npn4*4+1)) 
 
      CALL fact
      CALL signum
      lmax=2*nmax
      l1max=lmax+1
      ci=(0d0,1d0)
      cim(1)=ci
      DO 2 i=2,nmax
         cim(i)=cim(i-1)*ci
    2 CONTINUE
      ssi(1)=1d0
      DO 3 i=1,lmax
         i1=i+1
         si=dfloat(2*i+1)
         ssi(i1)=si
         IF(i.LE.nmax) ssj(i)=dsqrt(si)
    3 CONTINUE
      ci=-ci
      DO 5 i=1,nmax
         si=ssj(i)
         cci=cim(i)
         DO 4 j=1,nmax
            sj=1d0/ssj(j)
            ccj=cim(j)*sj/cci
            fr(j,i)=ccj
            fi(j,i)=ccj*ci
            ff(j,i)=si*sj
    4    CONTINUE
    5 CONTINUE
      nmax1=nmax+1
 
C *****  CALCULATION OF THE ARRAYS B1 AND B2  *****
 
      k1=1
      k2=0
      k3=0
      k4=1
      k5=1
      k6=2
 
C     PRINT 3300, B1,B2
 3300 FORMAT (' B1 AND B2')
      DO 100 n=1,nmax
 
C *****  CALCULATION OF THE ARRAYS T1 AND T2  *****
 
 
         DO 10 nn=1,nmax
            m1max=min0(n,nn)+1
            DO 6 m1=1,m1max
               m=m1-1
               l1=npn6+m
               tt1=tr11(m1,n,nn)
               tt2=tr12(m1,n,nn)
               tt3=tr21(m1,n,nn)
               tt4=tr22(m1,n,nn)
               tt5=ti11(m1,n,nn)
               tt6=ti12(m1,n,nn)
               tt7=ti21(m1,n,nn)
               tt8=ti22(m1,n,nn)
               t1=tt1+tt2
               t2=tt3+tt4
               t3=tt5+tt6
               t4=tt7+tt8
               tr1(l1,nn)=t1+t2
               tr2(l1,nn)=t1-t2
               ti1(l1,nn)=t3+t4
               ti2(l1,nn)=t3-t4
               IF(m.EQ.0) GO TO 6
               l1=npn6-m
               t1=tt1-tt2
               t2=tt3-tt4
               t3=tt5-tt6
               t4=tt7-tt8
               tr1(l1,nn)=t1-t2
               tr2(l1,nn)=t1+t2
               ti1(l1,nn)=t3-t4
               ti2(l1,nn)=t3+t4
    6       CONTINUE
   10    CONTINUE
 
C  *****  END OF THE CALCULATION OF THE ARRAYS T1 AND T2  *****
 
         nn1max=nmax1+n
         DO 40 nn1=1,nn1max
            n1=nn1-1
 
C  *****  CALCULATION OF THE ARRAYS A1 AND A2  *****
 
            CALL ccg(n,n1,nmax,k1,k2,g1)
            nnmax=min0(nmax,n1+n)
            nnmin=max0(1,iabs(n-n1))
            kn=n+nn1
            DO 15 nn=nnmin,nnmax
               nnn=nn+1
               sig=ssign(kn+nn)
               m1max=min0(n,nn)+npn6
               aar1=0d0
               aar2=0d0
               aai1=0d0
               aai2=0d0
               DO 13 m1=npn6,m1max
                  m=m1-npn6
                  sss=g1(m1,nnn)
                  rr1=tr1(m1,nn)
                  ri1=ti1(m1,nn)
                  rr2=tr2(m1,nn)
                  ri2=ti2(m1,nn)
                  IF(m.EQ.0) GO TO 12
                  m2=npn6-m
                  rr1=rr1+tr1(m2,nn)*sig
                  ri1=ri1+ti1(m2,nn)*sig
                  rr2=rr2+tr2(m2,nn)*sig
                  ri2=ri2+ti2(m2,nn)*sig
   12             aar1=aar1+sss*rr1
                  aai1=aai1+sss*ri1
                  aar2=aar2+sss*rr2
                  aai2=aai2+sss*ri2
   13          CONTINUE
               xr=fr(nn,n)
               xi=fi(nn,n)
               ar1(nn)=aar1*xr-aai1*xi
               ai1(nn)=aar1*xi+aai1*xr
               ar2(nn)=aar2*xr-aai2*xi
               ai2(nn)=aar2*xi+aai2*xr
   15       CONTINUE
 
C  *****  END OF THE CALCULATION OF THE ARRAYS A1 AND A2 ****
 
            CALL ccg(n,n1,nmax,k3,k4,g2)
            m1=max0(-n1+1,-n)
            m2=min0(n1+1,n)
            m1max=m2+npn6
            m1min=m1+npn6
            DO 30 m1=m1min,m1max
               bbr1=0d0
               bbi1=0d0
               bbr2=0d0
               bbi2=0d0
               DO 25 nn=nnmin,nnmax
                  nnn=nn+1
                  sss=g2(m1,nnn)
                  bbr1=bbr1+sss*ar1(nn)
                  bbi1=bbi1+sss*ai1(nn)
                  bbr2=bbr2+sss*ar2(nn)
                  bbi2=bbi2+sss*ai2(nn)
   25          CONTINUE
               b1r(nn1,m1,n)=bbr1
               b1i(nn1,m1,n)=bbi1
               b2r(nn1,m1,n)=bbr2
               b2i(nn1,m1,n)=bbi2
   30       CONTINUE
   40    CONTINUE
  100 CONTINUE
 
C  *****  END OF THE CALCULATION OF THE ARRAYS B1 AND B2 ****
 
C  *****  CALCULATION OF THE ARRAYS D1,D2,D3,D4, AND D5  *****
 
C     PRINT 3301
 3301 FORMAT(' D1, D2, ...')
      DO 200 n=1,nmax
         DO 190 nn=1,nmax
            m1=min0(n,nn)
            m1max=npn6+m1
            m1min=npn6-m1
            nn1max=nmax1+min0(n,nn)
            DO 180 m1=m1min,m1max
               m=m1-npn6
               nn1min=iabs(m-1)+1
               dd1=0d0
               dd2=0d0
               DO 150 nn1=nn1min,nn1max
                  xx=ssi(nn1)
                  x1=b1r(nn1,m1,n)
                  x2=b1i(nn1,m1,n)
                  x3=b1r(nn1,m1,nn)
                  x4=b1i(nn1,m1,nn)
                  x5=b2r(nn1,m1,n)
                  x6=b2i(nn1,m1,n)
                  x7=b2r(nn1,m1,nn)
                  x8=b2i(nn1,m1,nn)
                  dd1=dd1+xx*(x1*x3+x2*x4)
                  dd2=dd2+xx*(x5*x7+x6*x8)
  150          CONTINUE
               d1(m1,nn,n)=dd1
               d2(m1,nn,n)=dd2
  180       CONTINUE
            mmax=min0(n,nn+2)
            mmin=max0(-n,-nn+2)
            m1max=npn6+mmax
            m1min=npn6+mmin
            DO 186 m1=m1min,m1max
               m=m1-npn6
               nn1min=iabs(m-1)+1
               dd3=0d0
               dd4=0d0
               dd5r=0d0
               dd5i=0d0
               m2=-m+2+npn6
               DO 183 nn1=nn1min,nn1max
                  xx=ssi(nn1)
                  x1=b1r(nn1,m1,n)
                  x2=b1i(nn1,m1,n)
                  x3=b2r(nn1,m1,n)
                  x4=b2i(nn1,m1,n)
                  x5=b1r(nn1,m2,nn)
                  x6=b1i(nn1,m2,nn)
                  x7=b2r(nn1,m2,nn)
                  x8=b2i(nn1,m2,nn)
                  dd3=dd3+xx*(x1*x5+x2*x6)
                  dd4=dd4+xx*(x3*x7+x4*x8)
                  dd5r=dd5r+xx*(x3*x5+x4*x6)
                  dd5i=dd5i+xx*(x4*x5-x3*x6)
  183          CONTINUE
               d3(m1,nn,n)=dd3
               d4(m1,nn,n)=dd4
               d5r(m1,nn,n)=dd5r
               d5i(m1,nn,n)=dd5i
  186       CONTINUE
  190    CONTINUE
  200 CONTINUE
 
C  *****  END OF THE CALCULATION OF THE D-ARRAYS *****
 
C  *****  CALCULATION OF THE EXPANSION COEFFICIENTS *****
 
C     PRINT 3303
 3303 FORMAT (' G1, G2, ...')
 
      dk=lam*lam/(4d0*csca*dacos(-1d0))
      DO 300 l1=1,l1max
         g1l=0d0
         g2l=0d0
         g3l=0d0
         g4l=0d0
         g5lr=0d0
         g5li=0d0
         l=l1-1
         sl=ssi(l1)*dk
         DO 290 n=1,nmax
            nnmin=max0(1,iabs(n-l))
            nnmax=min0(nmax,n+l)
            IF(nnmax.LT.nnmin) GO TO 290
            CALL ccg(n,l,nmax,k1,k2,g1)
            IF(l.GE.2) CALL ccg(n,l,nmax,k5,k6,g2)
            nl=n+l
            DO 280  nn=nnmin,nnmax
               nnn=nn+1
               mmax=min0(n,nn)
               m1min=npn6-mmax
               m1max=npn6+mmax
               si=ssign(nl+nnn)
               dm1=0d0
               dm2=0d0
               DO 270 m1=m1min,m1max
                  m=m1-npn6
                  IF(m.GE.0) sss1=g1(m1,nnn)
                  IF(m.LT.0) sss1=g1(npn6-m,nnn)*si
                  dm1=dm1+sss1*d1(m1,nn,n)
                  dm2=dm2+sss1*d2(m1,nn,n)
  270          CONTINUE
               ffn=ff(nn,n)
               sss=g1(npn6+1,nnn)*ffn
               g1l=g1l+sss*dm1
               g2l=g2l+sss*dm2*si
               IF(l.LT.2) GO TO 280
               dm3=0d0
               dm4=0d0
               dm5r=0d0
               dm5i=0d0
               mmax=min0(n,nn+2)
               mmin=max0(-n,-nn+2)
               m1max=npn6+mmax
               m1min=npn6+mmin
               DO 275 m1=m1min,m1max
                  m=m1-npn6
                  sss1=g2(npn6-m,nnn)
                  dm3=dm3+sss1*d3(m1,nn,n)
                  dm4=dm4+sss1*d4(m1,nn,n)
                  dm5r=dm5r+sss1*d5r(m1,nn,n)
                  dm5i=dm5i+sss1*d5i(m1,nn,n)
  275          CONTINUE
               g5lr=g5lr-sss*dm5r
               g5li=g5li-sss*dm5i
               sss=g2(npn4,nnn)*ffn
               g3l=g3l+sss*dm3
               g4l=g4l+sss*dm4*si
  280       CONTINUE
  290    CONTINUE
         g1l=g1l*sl
         g2l=g2l*sl
         g3l=g3l*sl
         g4l=g4l*sl
         g5lr=g5lr*sl
         g5li=g5li*sl
         alf1(l1)=g1l+g2l
         alf2(l1)=g3l+g4l
         alf3(l1)=g3l-g4l
         alf4(l1)=g1l-g2l
         bet1(l1)=g5lr*2d0
         bet2(l1)=g5li*2d0
         lmax=l
         IF(abs(g1l).LT.1d-6) GO TO 500
  300 CONTINUE
  500 CONTINUE
      RETURN
      END
 
C****************************************************************
 
C   CALCULATION OF THE QUANTITIES F(N+1)=0.5*LN(N!)
C   0.LE.N.LE.899
 
      SUBROUTINE fact
      real*8 f(900)
      COMMON /fac/ f
      f(1)=0d0
      f(2)=0d0
      DO 2 i=3,900
         i1=i-1
         f(i)=f(i1)+0.5d0*dlog(dfloat(i1))
    2 CONTINUE
      RETURN
      END
 
C************************************************************
 
C   CALCULATION OF THE ARRAY SSIGN(N+1)=SIGN(N)
C   0.LE.N.LE.899
 
      SUBROUTINE signum
      real*8 ssign(900)
      COMMON /ss/ ssign
      ssign(1)=1d0
      DO 2 n=2,899 
         ssign(n)=-ssign(n-1)
    2 CONTINUE
      RETURN
      END
 
C******************************************************************
C
C   CALCULATION OF CLEBSCH-GORDAN COEFFICIENTS
C   (N,M:N1,M1/NN,MM)
C   FOR GIVEN N AND N1. M1=MM-M, INDEX MM IS FOUND FROM M AS
C   MM=M*K1+K2
C
C   INPUT PARAMETERS :  N,N1,NMAX,K1,K2
C                               N.LE.NMAX
C                               N.GE.1
C                               N1.GE.0
C                               N1.LE.N+NMAX
C   OUTPUT PARAMETERS : GG(M+NPN6,NN+1) - ARRAY OF THE CORRESPONDING
C                                       COEFFICIENTS
C                               /M/.LE.N
C                               /M1/=/M*(K1-1)+K2/.LE.N1
C                               NN.LE.MIN(N+N1,NMAX)
C                               NN.GE.MAX(/MM/,/N-N1/)
C   IF K1=1 AND K2=0, THEN 0.LE.M.LE.N
 
 
      SUBROUTINE ccg(N,N1,NMAX,K1,K2,GG)
      include 'tmd.par.f'
      IMPLICIT REAL*8 (a-h,o-z)
      real*8 gg(npl1,npn6),cd(0:npn5),cu(0:npn5)
      IF(nmax.LE.npn4.
     &   and.0.LE.n1.
     &   and.n1.LE.nmax+n.
     &   and.n.GE.1.
     &   and.n.LE.nmax) GO TO 1
      print 5001
      stop
 5001 FORMAT(' ERROR IN SUBROUTINE CCG')
    1 nnf=min0(n+n1,nmax)
      min=npn6-n
      mf=npn6+n
      IF(k1.EQ.1.AND.k2.EQ.0) min=npn6
      DO 100 mind=min,mf
         m=mind-npn6
         mm=m*k1+k2
         m1=mm-m
         IF(iabs(m1).GT.n1) GO TO 90
         nnl=max0(iabs(mm),iabs(n-n1))
         IF(nnl.GT.nnf) GO TO 90
         nnu=n+n1
         nnm=(nnu+nnl)*0.5d0
         IF (nnu.EQ.nnl) nnm=nnl
         CALL ccgin(n,n1,m,mm,c)
         cu(nnl)=c  
         IF (nnl.EQ.nnf) GO TO 50
         c2=0d0
         c1=c
         DO 7 nn=nnl+1,min0(nnm,nnf)
            a=dfloat((nn+mm)*(nn-mm)*(n1-n+nn))
            a=a*dfloat((n-n1+nn)*(n+n1-nn+1)*(n+n1+nn+1))
            a=dfloat(4*nn*nn)/a
            a=a*dfloat((2*nn+1)*(2*nn-1))
            a=dsqrt(a)
            b=0.5d0*dfloat(m-m1)
            d=0d0
            IF(nn.EQ.1) GO TO 5
            b=dfloat(2*nn*(nn-1))
            b=dfloat((2*m-mm)*nn*(nn-1)-mm*n*(n+1)+
     &               mm*n1*(n1+1))/b
            d=dfloat(4*(nn-1)*(nn-1))
            d=d*dfloat((2*nn-3)*(2*nn-1))
            d=dfloat((nn-mm-1)*(nn+mm-1)*(n1-n+nn-1))/d
            d=d*dfloat((n-n1+nn-1)*(n+n1-nn+2)*(n+n1+nn))
            d=dsqrt(d)
    5       c=a*(b*c1-d*c2)
            c2=c1
            c1=c
            cu(nn)=c
    7    CONTINUE
         IF (nnf.LE.nnm) GO TO 50
         CALL direct(n,m,n1,m1,nnu,mm,c)
         cd(nnu)=c
         IF (nnu.EQ.nnm+1) GO TO 50
         c2=0d0
         c1=c
         DO 12 nn=nnu-1,nnm+1,-1
            a=dfloat((nn-mm+1)*(nn+mm+1)*(n1-n+nn+1))
            a=a*dfloat((n-n1+nn+1)*(n+n1-nn)*(n+n1+nn+2))
            a=dfloat(4*(nn+1)*(nn+1))/a
            a=a*dfloat((2*nn+1)*(2*nn+3))
            a=dsqrt(a)
            b=dfloat(2*(nn+2)*(nn+1))
            b=dfloat((2*m-mm)*(nn+2)*(nn+1)-mm*n*(n+1)
     &               +mm*n1*(n1+1))/b
            d=dfloat(4*(nn+2)*(nn+2))
            d=d*dfloat((2*nn+5)*(2*nn+3))
            d=dfloat((nn+mm+2)*(nn-mm+2)*(n1-n+nn+2))/d
            d=d*dfloat((n-n1+nn+2)*(n+n1-nn-1)*(n+n1+nn+3))
            d=dsqrt(d)
            c=a*(b*c1-d*c2)
            c2=c1
            c1=c
            cd(nn)=c
   12    CONTINUE
   50    DO 9 nn=nnl,nnf
            IF (nn.LE.nnm) gg(mind,nn+1)=cu(nn)
            IF (nn.GT.nnm) gg(mind,nn+1)=cd(nn)
c           WRITE (6,*) N,M,N1,M1,NN,MM,GG(MIND,NN+1)
    9    CONTINUE
   90    CONTINUE
  100 CONTINUE
      RETURN
      END
 
C*********************************************************************
 
      SUBROUTINE direct (N,M,N1,M1,NN,MM,C)
      IMPLICIT REAL*8 (a-h,o-z)
      real*8 f(900)
      COMMON /fac/ f
      c=f(2*n+1)+f(2*n1+1)+f(n+n1+m+m1+1)+f(n+n1-m-m1+1)    
      c=c-f(2*(n+n1)+1)-f(n+m+1)-f(n-m+1)-f(n1+m1+1)-f(n1-m1+1)
      c=dexp(c)
      RETURN
      END
 
C*********************************************************************
C
C   CALCULATION OF THE CLEBCSH-GORDAN COEFFICIENTS
C   G=(N,M:N1,MM-M/NN,MM)
C   FOR GIVEN N,N1,M,MM, WHERE NN=MAX(/MM/,/N-N1/)
C                               /M/.LE.N
C                               /MM-M/.LE.N1
C                               /MM/.LE.N+N1
 
      SUBROUTINE ccgin(N,N1,M,MM,G)
      IMPLICIT REAL*8 (a-h,o-z)
      real*8 f(900),ssign(900)
      COMMON /ss/ ssign
      COMMON /fac/ f
      m1=mm-m
      IF(n.GE.iabs(m).
     &   and.n1.GE.iabs(m1).
     &   and.iabs(mm).LE.(n+n1)) GO TO 1
      print 5001
      stop
 5001 FORMAT(' ERROR IN SUBROUTINE CCGIN')
    1 IF (iabs(mm).GT.iabs(n-n1)) GO TO 100
      l1=n
      l2=n1
      l3=m
      IF(n1.LE.n) GO TO 50
      k=n
      n=n1
      n1=k
      k=m
      m=m1
      m1=k
   50 n2=n*2
      m2=m*2
      n12=n1*2
      m12=m1*2
      g=ssign(n1+m1+1)
     & *dexp(f(n+m+1)+f(n-m+1)+f(n12+1)+f(n2-n12+2)-f(n2+2)
     &       -f(n1+m1+1)-f(n1-m1+1)-f(n-n1+mm+1)-f(n-n1-mm+1))
      n=l1
      n1=l2
      m=l3
      RETURN
  100 a=1d0
      l1=m
      l2=mm
      IF(mm.GE.0) GO TO 150
      mm=-mm
      m=-m
      m1=-m1
      a=ssign(mm+n+n1+1)
  150 g=a*ssign(n+m+1)
     &   *dexp(f(2*mm+2)+f(n+n1-mm+1)+f(n+m+1)+f(n1+m1+1)
     &        -f(n+n1+mm+2)-f(n-n1+mm+1)-f(-n+n1+mm+1)-f(n-m+1)
     &        -f(n1-m1+1))
      m=l1
      mm=l2
      RETURN
      END
 
C********************************************************************
 
c********************************************************************
 
C********************************************************************
 
C********************************************************************
 
C  COMPUTATION OF R1 AND R2 FOR A POWER LAW SIZE DISTRIBUTION WITH
C  EFFECTIVE RADIUS A AND EFFECTIVE VARIANCE B
 
      SUBROUTINE power (A,B,R1,R2)
      IMPLICIT REAL*8 (a-h,o-z)
      EXTERNAL f
      COMMON aa,bb
      aa=a
      bb=b
      ax=1d-5
      bx=a-1d-5
      r1=zeroin(ax,bx,f,0d0)
      r2=(1d0+b)*2d0*a-r1
      RETURN
      END
 
C***********************************************************************
 
      DOUBLE PRECISION FUNCTION zeroin (AX,BX,F,TOL)
      IMPLICIT REAL*8 (a-h,o-z)
      eps=1d0
   10 eps=0.5d0*eps
      tol1=1d0+eps
      IF (tol1.GT.1d0) GO TO 10
   15 a=ax
      b=bx
      fa=f(a)
      fb=f(b)
   20 c=a
      fc=fa
      d=b-a
      e=d
   30 IF (abs(fc).GE.abs(fb)) GO TO 40
   35 a=b
      b=c
      c=a
      fa=fb
      fb=fc
      fc=fa
   40 tol1=2d0*eps*abs(b)+0.5d0*tol
      xm=0.5d0*(c-b)
      IF (abs(xm).LE.tol1) GO TO 90
   44 IF (fb.EQ.0d0) GO TO 90
   45 IF (abs(e).LT.tol1) GO TO 70
   46 IF (abs(fa).LE.abs(fb)) GO TO 70
   47 IF (a.NE.c) GO TO 50
   48 s=fb/fa
      p=2d0*xm*s
      q=1d0-s
      GO TO 60
   50 q=fa/fc
      r=fb/fc
      s=fb/fa
      p=s*(2d0*xm*q*(q-r)-(b-a)*(r-1d0))
      q=(q-1d0)*(r-1d0)*(s-1d0)
   60 IF (p.GT.0d0) q=-q
      p=abs(p)
      IF ((2d0*p).GE.(3d0*xm*q-abs(tol1*q))) GO TO 70
   64 IF (p.GE.abs(0.5d0*e*q)) GO TO 70
   65 e=d
      d=p/q
      GO TO 80
   70 d=xm
      e=d
   80 a=b
      fa=fb
      IF (abs(d).GT.tol1) b=b+d
      IF (abs(d).LE.tol1) b=b+dsign(tol1,xm)
      fb=f(b)
      IF ((fb*(fc/abs(fc))).GT.0d0) GO TO 20
   85 GO TO 30
   90 zeroin=b
      RETURN
      END
 
C***********************************************************************
 
      DOUBLE PRECISION FUNCTION f(R1)
      IMPLICIT REAL*8 (a-h,o-z)
      COMMON a,b
      r2=(1d0+b)*2d0*a-r1
      f=(r2-r1)/dlog(r2/r1)-a
      RETURN
      END
 
C**********************************************************************
C    CALCULATION OF POINTS AND WEIGHTS OF GAUSSIAN QUADRATURE         *
C    FORMULA. IF IND1 = 0 - ON INTERVAL (-1,1), IF IND1 = 1 - ON      *
C    INTERVAL  (0,1). IF  IND2 = 1 RESULTS ARE PRINTED.               *
C    N - NUMBER OF POINTS                                             *
C    Z - DIVISION POINTS                                              *
C    W - WEIGHTS                                                      *
C**********************************************************************
 
C****************************************************************
 
      SUBROUTINE distrb (NNK,YY,WY,NDISTR,AA,BB,GAM,R1,R2,REFF,                 
     &                   VEFF,PI)                                               
      IMPLICIT REAL*8 (a-h,o-z)                                                 
      real*8 yy(nnk),wy(nnk)                                                    
      IF (ndistr.EQ.2) GO TO 100                                                
      IF (ndistr.EQ.3) GO TO 200                                                
      IF (ndistr.EQ.4) GO TO 300                                                
      IF (ndistr.EQ.5) GO TO 360
      print 1001,aa,bb,gam                                                      
 1001 FORMAT('MODIFIED GAMMA DISTRIBUTION, alpha=',f6.4,'  r_c=',               
     &  f6.4,'  gamma=',f6.4)                                                   
      a2=aa/gam                                                                 
      db=1d0/bb
      DO 50 i=1,nnk                                                             
         x=yy(i)                                                             
         y=x**aa                                                                
         x=x*db
         y=y*dexp(-a2*(x**gam))                                                 
         wy(i)=wy(i)*y                                                       
   50 CONTINUE                                                                  
      GO TO 400                                                                 
  100 print 1002,aa,bb                                                          
 1002 FORMAT('LOG-NORMAL DISTRIBUTION, r_g=',f8.4,                
     &       '  [ln(sigma_g)]**2=', f6.4)          
      da=1d0/aa                                                                 
      DO 150 i=1,nnk                                                            
         x=yy(i)                                                                
         y=dlog(x*da)                                                          
         y=dexp(-y*y*0.5d0/bb)/x                                             
         wy(i)=wy(i)*y                                                          
  150 CONTINUE                                                                  
      GO TO 400                                                                 
  200 print 1003                                                                
 1003 FORMAT('POWER LAW DISTRIBUTION OF HANSEN & TRAVIS 1974')                 
      DO 250 i=1,nnk                                                            
         x=yy(i)                                                                
         wy(i)=wy(i)/(x*x*x)                                                 
  250 CONTINUE                                                                  
      GO TO 400                                                                 
  300 print 1004,aa,bb                                                          
 1004 FORMAT ('GAMMA DISTRIBUTION,  a=',f6.3,'  b=',f6.4)
      b2=(1d0-3d0*bb)/bb                                                        
      dab=1d0/(aa*bb)                                                          
      DO 350 i=1,nnk                                                            
         x=yy(i)                                                                
         x=(x**b2)*dexp(-x*dab)                                                 
         wy(i)=wy(i)*x                                                       
  350 CONTINUE                                                                  
      GO TO 400                                                                 
  360 print 1005,bb
 1005 FORMAT ('MODIFIED POWER LAW DISTRIBUTION,  alpha=',d10.4)
      DO 370 i=1,nnk
         x=yy(i)
         IF (x.LE.r1) wy(i)=wy(i)
         IF (x.GT.r1) wy(i)=wy(i)*(x/r1)**bb
  370 CONTINUE
  400 CONTINUE                                                                  
      sum=0d0
      DO 450 i=1,nnk
         sum=sum+wy(i)
  450 CONTINUE
      sum=1d0/sum
      DO 500 i=1,nnk
         wy(i)=wy(i)*sum
  500 CONTINUE
      g=0d0
      DO 550 i=1,nnk
         x=yy(i)
         g=g+x*x*wy(i)
  550 CONTINUE
      reff=0d0
      DO 600 i=1,nnk
         x=yy(i)
         reff=reff+x*x*x*wy(i)
  600 CONTINUE
      reff=reff/g
      veff=0d0
      DO 650 i=1,nnk
         x=yy(i)
         xi=x-reff
         veff=veff+xi*xi*x*x*wy(i)
  650 CONTINUE
      veff=veff/(g*reff*reff)
      RETURN                                                                    
      END                                                                       
 
C*************************************************************
 
      SUBROUTINE hovenr(L1,A1,A2,A3,A4,B1,B2)
      IMPLICIT REAL*8 (a-h,o-z)
      real*8 a1(l1),a2(l1),a3(l1),a4(l1),b1(l1),b2(l1)
      DO 100 l=1,l1
         kontr=1
         ll=l-1
         dl=dfloat(ll)*2d0+1d0
         ddl=dl*0.48d0
         aa1=a1(l)
         aa2=a2(l)
         aa3=a3(l)
         aa4=a4(l)
         bb1=b1(l)
         bb2=b2(l)
         IF(ll.GE.1.AND.abs(aa1).GE.dl) kontr=2
         IF(abs(aa2).GE.dl) kontr=2
         IF(abs(aa3).GE.dl) kontr=2
         IF(abs(aa4).GE.dl) kontr=2
         IF(abs(bb1).GE.ddl) kontr=2
         IF(abs(bb2).GE.ddl) kontr=2
         IF(kontr.EQ.2) print 3000,ll
         c=-0.1d0
         DO 50 i=1,11
            c=c+0.1d0
            cc=c*c
            c1=cc*bb2*bb2
            c2=c*aa4
            c3=c*aa3
            IF((dl-c*aa1)*(dl-c*aa2)-cc*bb1*bb1.LE.-1d-4) kontr=2
            IF((dl-c2)*(dl-c3)+c1.LE.-1d-4) kontr=2
            IF((dl+c2)*(dl-c3)-c1.LE.-1d-4) kontr=2
            IF((dl-c2)*(dl+c3)-c1.LE.-1d-4) kontr=2
            IF(kontr.EQ.2) print 4000,ll,c
   50    CONTINUE
  100 CONTINUE
      IF(kontr.EQ.1) print 2000
 2000 FORMAT('TEST OF VAN DER MEE & HOVENIER IS SATISFIED')
 3000 FORMAT('TEST OF VAN DER MEE & HOVENIER IS NOT SATISFIED, L=',i3)
 4000 FORMAT('TEST OF VAN DER MEE & HOVENIER IS NOT SATISFIED, L=',i3,
     & '   A=',d9.2)
      RETURN
      END
 
C****************************************************************
 
C    CALCULATION OF THE SCATTERING MATRIX FOR GIVEN EXPANSION
C    COEFFICIENTS
 
C    A1,...,B2 - EXPANSION COEFFICIENTS
C    LMAX - NUMBER OF COEFFICIENTS MINUS 1
C    N - NUMBER OF SCATTERING ANGLES
C        THE CORRESPONDING SCATTERING ANGLES ARE GIVEN BY
C        180*(I-1)/(N-1) (DEGREES), WHERE I NUMBERS THE ANGLES
 
      SUBROUTINE matr(A1,A2,A3,A4,B1,B2,LMAX,NPNA,F)
      include 'tmd.par.f'
      IMPLICIT REAL*8 (a-h,o-z)
      real*8 a1(npl),a2(npl),a3(npl),a4(npl),b1(npl),b2(npl),f(npna,4,4)
      n=npna
      dn=1d0/dfloat(n-1)
      da=dacos(-1d0)*dn
      db=180d0*dn
      l1max=lmax+1
      print 1000
 1000 FORMAT(' ')
      print 1001
 1001 FORMAT(' ',2x,'S',6x,'ALPHA1',6x,'ALPHA2',6x,'ALPHA3',
     &       6x,'ALPHA4',7x,'BETA1',7x,'BETA2')
      DO 10 l1=1,l1max
         l=l1-1
         print 1002,l,a1(l1),a2(l1),a3(l1),a4(l1),b1(l1),b2(l1)
   10 CONTINUE
 1002 FORMAT(' ',i3,6f12.5)
      tb=-db
      taa=-da
      print 1000
      print 1003
 1003 FORMAT(' ',5x,'<',8x,'F11',8x,'F22',8x,'F33',
     & 8x,'F44',8x,'F12',8x,'F34')
      d6=dsqrt(6d0)*0.25d0
      DO 500 i1=1,n
         taa=taa+da
         tb=tb+db
         u=dcos(taa)
         f11=0d0
         f2=0d0
         f3=0d0
         f44=0d0
         f12=0d0
         f34=0d0
         p1=0d0
         p2=0d0
         p3=0d0
         p4=0d0
         pp1=1d0
         pp2=0.25d0*(1d0+u)*(1d0+u)
         pp3=0.25d0*(1d0-u)*(1d0-u)
         pp4=d6*(u*u-1d0)
         DO 400 l1=1,l1max
            l=l1-1
            dl=dfloat(l)
            dl1=dfloat(l1)
            f11=f11+a1(l1)*pp1
            f44=f44+a4(l1)*pp1
            IF(l.EQ.lmax) GO TO 350
            pl1=dfloat(2*l+1)
            p=(pl1*u*pp1-dl*p1)/dl1
            p1=pp1
            pp1=p
  350       IF(l.LT.2) GO TO 400
            f2=f2+(a2(l1)+a3(l1))*pp2
            f3=f3+(a2(l1)-a3(l1))*pp3
            f12=f12+b1(l1)*pp4
            f34=f34+b2(l1)*pp4
            IF(l.EQ.lmax) GO TO 400
            pl2=dfloat(l*l1)*u
            pl3=dfloat(l1*(l*l-4))
            pl4=1d0/dfloat(l*(l1*l1-4))
            p=(pl1*(pl2-4d0)*pp2-pl3*p2)*pl4
            p2=pp2
            pp2=p
            p=(pl1*(pl2+4d0)*pp3-pl3*p3)*pl4
            p3=pp3
            pp3=p
            p=(pl1*u*pp4-dsqrt(dfloat(l*l-4))*p4)/dsqrt(dfloat(l1*l1-4))
            p4=pp4
            pp4=p
  400    CONTINUE
         f22=(f2+f3)*0.5d0
         f33=(f2-f3)*0.5d0
C        F22=F22/F11
C        F33=F33/F11
C        F44=F44/F11
C        F12=-F12/F11
C        F34=F34/F11
         print 1004,tb,f11,f22,f33,f44,f12,f34
         f(i1,1,1)=f11
         f(i1,1,2)=f12
         f(i1,1,3)=0
         f(i1,1,4)=0
         f(i1,2,1)=f12
         f(i1,2,2)=f22
         f(i1,2,3)=0
         f(i1,2,4)=0
         f(i1,3,1)=0
         f(i1,3,2)=0
         f(i1,3,3)=f33
         f(i1,3,4)=f34
         f(i1,4,1)=0
         f(i1,4,2)=0
         f(i1,4,3)=-f34
         f(i1,4,4)=f44
  500 CONTINUE
      print 1000 
 1004 FORMAT(' ',f6.2,6f11.4)
      RETURN
      END