OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
tm_tmd.lp.f
Go to the documentation of this file.
1 C This is NOT the original code by the authors listed below!
2 C Modified 2009-2013 by Jussi Leinonen (jsleinonen@gmail.com)
3 C to be compatible with the Python extension interface.
4 C The requirement of non-profit use only has been dropped from
5 C this code with the permission of M. I. Mishchenko; thus the
6 C license has been made open source compatible.
7 
8 C New release including the LAPACK matrix inversion procedure.
9 C We thank Cory Davis (University of Edinburgh) for pointing
10 C out the possibility of replacing the proprietary NAG matrix
11 C inversion routine by the public-domain LAPACK equivalent.
12 
13 C CALCULATION OF LIGHT SCATTERING BY POLYDISPERSE, RANDOMLY
14 C ORIENTED PARTICLES OF IDENTICAL AXIALLY SYMMETRIC SHAPE
15 
16 C This version of the code uses DOUBLE PRECISION variables
17 C and must be used along with the accompanying files tm_tmd.par.f
18 C and lpd.f.
19 
20 C Last update 08/06/2005
21 
22 C The code has been developed by Michael Mishchenko at the NASA
23 C Goddard Institute for Space Studies, New York. This research
24 C was funded by the NASA Radiation Sciences Program.
25 
26 C The code can be used without limitations in any not-for-
27 C profit scientific research. We only request that in any
28 C publication using the code the source of the code be acknowledged
29 C and relevant references (see below) be made.
30 
31 C This version of the code is applicable to spheroids,
32 C Chebyshev particles, and finite circular cylinders.
33 
34 C The computational method is based on the Watermsn's T-matrix
35 C approach and is described in detail in the following papers:
36 C
37 C 1. M. I. Mishchenko, Light scattering by randomly oriented
38 C axially symmetric particles, J. Opt. Soc. Am. A,
39 C vol. 8, 871-882 (1991).
40 C
41 C 2. M. I. Mishchenko, Light scattering by size-shape
42 C distributions of randomly oriented axially symmetric
43 C particles of a size comparable to a wavelength,
44 C Appl. Opt., vol. 32, 4652-4666 (1993).
45 C
46 C 3. M. I. Mishchenko and L. D. Travis, T-matrix computations
47 C of light scattering by large spheroidal particles,
48 C Opt. Commun., vol. 109, 16-21 (1994).
49 C
50 C 4. M. I. Mishchenko, L. D. Travis, and A. Macke, Scattering
51 C of light by polydisperse, randomly oriented, finite
52 C circular cylinders, Appl. Opt., vol. 35, 4927-4940 (1996).
53 C
54 C 5. D. J. Wielaard, M. I. Mishchenko, A. Macke, and B. E. Carlson,
55 C Improved T-matrix computations for large, nonabsorbing and
56 C weakly absorbing nonspherical particles and comparison
57 C with geometrical optics approximation, Appl. Opt., vol. 36,
58 C 4305-4313 (1997).
59 C
60 C A general review of the T-matrix approach can be found in
61 C
62 C 6. M. I. Mishchenko, L. D. Travis, and D. W. Mackowski,
63 C T-matrix computations of light scattering by nonspherical
64 C particles: a review, J. Quant. Spectrosc. Radiat.
65 C Transfer, vol. 55, 535-575 (1996).
66 C
67 C The following paper provides a detailed user guide to the
68 C T-matrix code:
69 C
70 C 7. M. I. Mishchenko and L. D. Travis, Capabilities and
71 C limitations of a current FORTRAN implementation of the
72 C T-matrix method for randomly oriented, rotationally
73 C symmetric scatterers, J. Quant. Spectrosc. Radiat. Transfer,
74 C vol. 60, 309-324 (1998).
75 C
76 C These papers are available in the .pdf format at the web site
77 C
78 C http://www.giss.nasa.gov/~crmim/publications/
79 C
80 C or in hardcopy upon request from Michael Mishchenko
81 C Please e-mail your request to crmim@giss.nasa.gov.
82 C
83 C A comprehensive book "Scattering, Absorption, and Emission of
84 C Light by Small Particles" (Cambridge University Press, Cambridge,
85 C 2002) is also available in the .pdf format at the web site
86 C
87 C http://www.giss.nasa.gov/~crmim/books.html
88 
89 C Analytical averaging over particle orientations (Ref. 1) makes
90 C this method the fastest exact technique currently available.
91 C The use of an automatic convergence procedure
92 C (Ref. 2) makes the code convenient in massive computations.
93 C Ref. 4 describes features specific for finite cylinders as
94 C particles with sharp rectangular edges. Ref. 5 describes further
95 C numerical improvements.
96 
97 C The use of extended precision variables (Ref. 3) can
98 C significantly increase the maximal convergent equivalent-sphere
99 C size parameter and make it greater than 200 (depending on
100 C refractive index and aspect ratio). The extended-precision code
101 C is also available. However, the use of extended precision varibales
102 C results in a greater consumption of CPU time.
103 C On IBM RISC workstations, that code is approximately
104 C five times slower than this double-precision code. The
105 C CPU time difference between the double-precision and extended-
106 C precision codes can be larger on supercomputers.
107 
108 C This is the first part of the full T-matrix code. The second part,
109 C lpd.f, is completely independent of the firsti part. It contains no
110 C T-matrix-specific subroutines and can be compiled separately.
111 C The second part of the code replaces the previously implemented
112 C standard matrix inversion scheme based on Gaussian elimination
113 C by a scheme based on the LU factorization technique.
114 C As described in Ref. 5 above, the use of the LU factorization is
115 C especially beneficial for nonabsorbing or weakly absorbing particles.
116 C In this code we use the LAPACK implementation of the LU factorization
117 C scheme. LAPACK stands for Linear Algebra PACKage. The latter is
118 C publicly available at the following internet site:
119 C
120 C http://www.netlib.org/lapack/
121 
122 C INPUT PARAMETERS:
123 C
124 C RAT = 1 - particle size is specified in terms of the
125 C equal-volume-sphere radius
126 C RAT.NE.1 - particle size is specified in terms of the
127 C equal-surface-area-sphere radius
128 C NDISTR specifies the distribution of equivalent-sphere radii
129 C NDISTR = 1 - modified gamma distribution
130 C [Eq. (40) of Ref. 7]
131 C AXI=alpha
132 C B=r_c
133 C GAM=gamma
134 C NDISTR = 2 - log-normal distribution
135 C [Eq. 41) of Ref. 7]
136 C AXI=r_g
137 C B=[ln(sigma_g)]**2
138 C NDISTR = 3 - power law distribution
139 C [Eq. (42) of Ref. 7]
140 C AXI=r_eff (effective radius)
141 C B=v_eff (effective variance)
142 C Parameters R1 and R2 (see below) are calculated
143 C automatically for given AXI and B
144 C NDISTR = 4 - gamma distribution
145 C [Eq. (39) of Ref. 7]
146 C AXI=a
147 C B=b
148 C NDISTR = 5 - modified power law distribution
149 C [Eq. (24) in M. I. Mishchenko et al.,
150 C Bidirectional reflectance of flat,
151 C optically thick particulate laters: an efficient radiative
152 C transfer solution and applications to snow and soil surfaces,
153 C J. Quant. Spectrosc. Radiat. Transfer, Vol. 63, 409-432 (1999)].
154 C B=alpha
155 C
156 C The code computes NPNAX size distributions of the same type
157 C and with the same values of B and GAM in one run.
158 C The parameter AXI varies from AXMAX to AXMAX/NPNAX in steps of
159 C AXMAX/NPNAX. To compute a single size distribution, use
160 C NPNAX=1 and AXMAX equal to AXI of this size distribution.
161 C
162 C R1 and R2 - minimum and maximum equivalent-sphere radii
163 C in the size distribution. They are calculated automatically
164 C for the power law distribution with given AXI and B
165 C but must be specified for other distributions
166 C after the lines
167 C
168 C DO 600 IAX=1,NPNAX
169 C AXI=AXMAX-DAX*DFLOAT(IAX-1)
170 C
171 C in the main program.
172 C For the modified power law distribution (NDISTR=5), the
173 C minimum radius is 0, R2 is the maximum radius,
174 C and R1 is the intermediate radius at which the
175 C n(r)=const dependence is replaced by the power law
176 C dependence.
177 C
178 C NKMAX.LE.988 is such that NKMAX+2 is the
179 C number of Gaussian quadrature points used in
180 C integrating over the size distribution for particles with
181 C AXI=AXMAX. For particles with AXI=AXMAX-AXMAX/NPNAX,
182 C AXMAX-2*AXMAX/NPNAX, etc. the number of Gaussian points
183 C linearly decreases.
184 C For the modified power law distribution, the number
185 C of integration points on the interval [0,R1] is also
186 C equal to NKMAX.
187 C
188 C LAM - wavelength of light
189 C MRR and MRI - real and imaginary parts of the refractive
190 C index (MRI.GE.0)
191 C EPS and NP - specify the shape of the particles.
192 C For spheroids NP=-1 and EPS is the ratio of the
193 C horizontal to rotational axes. EPS is larger than
194 C 1 for oblate spheroids and smaller than 1 for
195 C prolate spheroids.
196 C For cylinders NP=-2 and EPS is the ratio of the
197 C diameter to the length.
198 C For Chebyshev particles NP must be positive and
199 C is the degree of the Chebyshev polynomial, while
200 C EPS is the deformation parameter
201 C [Eq. (33) of Ref. 7].
202 C DDELT - accuracy of the computations
203 C NPNA - number of equidistant scattering angles (from 0
204 C to 180 deg) for which the scattering matrix is
205 C calculated.
206 C NDGS - parameter controlling the number of division points
207 C in computing integrals over the particle surface.
208 C For compact particles, the recommended value is 2.
209 C For highly aspherical particles larger values (3, 4,...)
210 C may be necessary to obtain convergence.
211 C The code does not check convergence over this parameter.
212 C Therefore, control comparisons of results obtained with
213 C different NDGS-values are recommended.
214 
215 
216 C OUTPUT PARAMETERS:
217 C
218 C REFF and VEFF - effective radius and effective variance of
219 C the size distribution as defined by Eqs. (43)-(45) of
220 C Ref. 7.
221 C CEXT - extinction cross section per particle
222 C CSCA - scattering cross section per particle
223 C W - single scattering albedo
224 C <cos> - asymmetry parameter of the phase function
225 C ALPHA1,...,BETA2 - coefficients appearing in the expansions
226 C of the elements of the scattering matrix in
227 C generalized spherical functions
228 C [Eqs. (11)-(16) of Ref. 7].
229 C F11,...,F44 - elements of the normalized scattering matrix [as
230 C defined by Eqs. (1)-(3) of Ref. 7] versus scattering angle
231 
232 C Note that LAM, r_c, r_g, r_eff, a, R1, and R2 must
233 C be given in the same units of length, and that
234 C the dimension of CEXT and CSCA is that of LAM squared (e.g., if
235 C LAM and AXI are given in microns, then CEXT and CSCA are
236 C calculated in square microns).
237 
238 C The physical correctness of the computed results is tested using
239 C the general inequalities derived by van der Mee and Hovenier,
240 C Astron. Astrophys., vol. 228, 559-568 (1990). Although
241 C the message that the test of van der Mee and Hovenier is satisfied
242 C does not guarantee that the results are absolutely correct,
243 C the message that the test is not satisfied can mean that something
244 C is wrong.
245 
246 C The convergence of the T-matrix method for particles with
247 C different sizes, refractive indices, and aspect ratios can be
248 C dramatically different. Usually, large sizes and large aspect
249 C ratios cause problems. The user of this code
250 C should first experiment with different input parameters in
251 C order to get an idea of the range of applicability of this
252 C technique. Sometimes decreasing the aspect ratio
253 C from 3 to 2 can increase the maximum convergent equivalent-
254 C sphere size parameter by a factor of several (Ref. 7).
255 C The CPU time required rapidly increases with increasing ratio
256 C radius/wavelength and/or with increasing particle asphericity.
257 C This should be taken into account in planning massive computations.
258 C Using an optimizing compiler on IBM RISC workstations saves
259 C about 70% of CPU time.
260 
261 C Execution can be automatically terminated if dimensions of certain
262 C arrays are not big enough or if the convergence procedure decides
263 C that the accuracy of double precision variables is insufficient
264 C to obtain a converged T-matrix solution for given particles.
265 C In all cases, a message appears explaining the cause of termination.
266 
267 C The message
268 C "WARNING: W IS GREATER THAN 1"
269 C means that the single-scattering albedo exceeds the maximum
270 C possible value 1. If W is greater than 1 by more than
271 C DDELT, this message can be an indication of numerical
272 C instability caused by extreme values of particle parameters.
273 
274 C The message "WARNING: NGAUSS=NPNG1" means that convergence over
275 C the parameter NG (see Ref. 2) cannot be obtained for the NPNG1
276 C value specified in the PARAMETER statement in the file tm_tmd.par.f.
277 C Often this is not a serious problem, especially for compact
278 C particles.
279 
280 C Larger and/or more aspherical particles may require larger
281 C values of the parameters NPN1, NPN4, and NPNG1 in the file
282 C tm_tmd.par.f. It is recommended to keep NPN1=NPN4+25 and
283 C NPNG1=3*NPN1. Note that the memory requirement increases
284 C as the third power of NPN4. If the memory of
285 C a computer is too small to accomodate the code in its current
286 C setting, the parameters NPN1, NPN4, and NPNG1 should be
287 C decreased. However, this will decrease the maximum size parameter
288 C that can be handled by the code.
289 
290 C In some cases any increases of NPN1 will not make the T-matrix
291 C computations convergent. This means that the particle is just
292 C too "bad" (extreme size parameter and/or extreme aspect ratio
293 C and/or extreme refractive index; see Ref. 7).
294 C The main program contains several PRINT statements which are
295 C currently commentd out. If uncommented, these statements will
296 C produce numbers which show the convergence rate and can be
297 C used to determine whether T-matrix computations for given particle
298 C parameters will converge at all.
299 
300 C Some of the common blocks are used to save memory rather than
301 C to transfer data. Therefore, if a compiler produces a warning
302 C message that the lengths of a common block are different in
303 C different subroutines, this is not a real problem.
304 
305 C The recommended value of DDELT is 0.001. For bigger values,
306 C false convergence can be obtained.
307 
308 C In computations for spheres use EPS=1.000001 instead of EPS=1.
309 C The use of EPS=1 can cause overflows in some rare cases.
310 
311 C To calculate a monodisperse particle, use the options
312 C NPNAX=1
313 C AXMAX=R
314 C B=1D-1
315 C NKMAX=-1
316 C NDISTR=4
317 C ...
318 C DO 600 IAX=1,NPNAX
319 C AXI=AXMAX-DAX*DFLOAT(IAX-1)
320 C R1=0.9999999*AXI
321 C R2=1.0000001*AXI
322 C ...
323 C where R is the equivalent-sphere radius.
324 
325 C It is recommended to use the power law rather than the
326 C gamma size distribution, because in this case convergent solution
327 C can be obtained for larger REFF and VEFF assuming the same
328 C maximal R2 (Mishchenko and Travis, Appl. Opt., vol. 33, 7206-7225,
329 C 1994).
330 
331 C For some compilers, DACOS must be raplaced by DARCOS and DASIN
332 C by DARSIN.
333 
334 C If many different size distributions are computed and the
335 C refractive index is fixed, then another approach can be more
336 C efficient than running this code many times. Specifically,
337 C scattering results should be computed for monodisperse particles
338 C with sizes ranging from essentially zero to some maximum value
339 C with a small step size (say, 0.02 microns). These results
340 C should be stored on disk and can be used along with spline
341 C interpolation to compute scattering by particles with intermediate
342 C sizes. Scattering patterns for monodisperse nonspherical
343 C particles in random orientation are (much) smoother than for
344 C monodisperse spheres, and spline interpolation usually gives good
345 C results. In this way, averaging over any size distribution is a
346 C matter of seconds. For more on size averaging, see Refs. 2 and 4.
347 
348 C We would highly appreciate informing me of any problems encountered
349 C with this code. Please send your message to the following
350 C e-mail address: CRMIM@GISS.NASA.GOV.
351 
352 C WHILE THE COMPUTER PROGRAM HAS BEEN TESTED FOR A VARIETY OF CASES,
353 C IT IS NOT INCONCEIVABLE THAT IT CONTAINS UNDETECTED ERRORS. ALSO,
354 C INPUT PARAMETERS CAN BE USED WHICH ARE OUTSIDE THE ENVELOPE OF
355 C VALUES FOR WHICH RESULTS ARE COMPUTED ACCURATELY. FOR THIS REASON,
356 C THE AUTHORS AND THEIR ORGANIZATION DISCLAIM ALL LIABILITY FOR
357 C ANY DAMAGES THAT MAY RESULT FROM THE USE OF THE PROGRAM.
358 
359  SUBROUTINE calcrand(RAT,LAM,MRR,MRI,EPS,NP,DDELT,NDGS,
360  & NPNAX,AXMAX,B,GAM,NDISTR,NKMAX,NPNA,NCOEFF,
361  & REFF,VEFF,CEXTIN,CSCAT,WALB,ASYMM, ALPHA, BETA, FM, SZ)
362  IMPLICIT REAL*8 (a-h,o-z)
363  include 'tm_tmd.par.f'
364  real*8 reff(npnax),veff(npnax),cextin(npnax)
365  real*8 cscat(npnax),walb(npnax),asymm(npnax)
366  real*8 alpha(4,npl,npnax),beta(2,npl,npnax),fm(6,npna,npnax),sz(npnax)
367 Cf2py intent(in) rat
368 Cf2py intent(in) lam
369 Cf2py intent(in) mrr
370 Cf2py intent(in) mri
371 Cf2py intent(in) eps
372 Cf2py intent(in) np
373 Cf2py intent(in) ddelt
374 Cf2py intent(in) ndgs
375 Cf2py intent(in) npnax
376 Cf2py intent(in) axmax
377 Cf2py intent(in) b
378 Cf2py intent(in) gam
379 Cf2py intent(in) ndistr
380 Cf2py intent(in) nkmax
381 Cf2py intent(in) npna
382 Cf2py intent(in) ncoeff
383 Cf2py intent(out) reff
384 Cf2py intent(out) veff
385 Cf2py intent(out) cextin
386 Cf2py intent(out) cscat
387 Cf2py intent(out) walb
388 Cf2py intent(out) asymm
389 Cf2py intent(out) alpha
390 Cf2py intent(out) beta
391 Cf2py intent(out) f
392 
393  real*8 lam,mrr,mri,x(npng2),w(npng2),s(npng2),ss(npng2),
394  * an(npn1),r(npng2),dr(npng2),
395  * ddr(npng2),drr(npng2),dri(npng2),ann(npn1,npn1)
396  real*8 xg(1000),wg(1000),tr1(npn2,npn2),ti1(npn2,npn2),
397  & alph1(npl),alph2(npl),alph3(npl),alph4(npl),bet1(npl),
398  & bet2(npl),xg1(2000),wg1(2000),
399  & al1(npl),al2(npl),al3(npl),al4(npl),be1(npl),be2(npl)
400  real*4
401  & rt11(npn6,npn4,npn4),rt12(npn6,npn4,npn4),
402  & rt21(npn6,npn4,npn4),rt22(npn6,npn4,npn4),
403  & it11(npn6,npn4,npn4),it12(npn6,npn4,npn4),
404  & it21(npn6,npn4,npn4),it22(npn6,npn4,npn4)
405 
406 
407  COMMON /ct/ tr1,ti1
408  COMMON /tmat/ rt11,rt12,rt21,rt22,it11,it12,it21,it22
409  p=dacos(-1d0)
410 
411  DO i=1,npn6
412  DO j=1,npn4
413  DO k=1,npn4
414  rt11(i,j,k) = 0.d0
415  rt12(i,j,k) = 0.d0
416  rt21(i,j,k) = 0.d0
417  rt22(i,j,k) = 0.d0
418  rt11(i,j,k) = 0.d0
419  rt12(i,j,k) = 0.d0
420  rt21(i,j,k) = 0.d0
421  rt22(i,j,k) = 0.d0
422  ENDDO
423  ENDDO
424  ENDDO
425 
426 C OPEN FILES *******************************************************
427 
428 C OPEN (6,FILE='test')
429 C OPEN (10,FILE='tmatr.write')
430 
431 C INPUT DATA ********************************************************
432 
433 C RAT=0.5 D0
434 C NDISTR=3
435 C AXMAX=1D0
436 C NPNAX=2
437 C B=0.1D0
438 C GAM=0.5D0
439 C NKMAX=5
440 C EPS=2D0
441 C NP=-1
442 C LAM=0.5D0
443 C MRR=1.53 d0
444 C MRI=0.008D0
445 C DDELT=0.001D0
446 C NPNA=19
447 C NDGS=2
448 
449  ncheck=0
450  IF (np.EQ.-1.OR.np.EQ.-2) ncheck=1
451  IF (np.GT.0.AND.(-1)**np.EQ.1) ncheck=1
452  WRITE (6,5454) ncheck
453  5454 FORMAT ('NCHECK=',i1)
454  dax=axmax/npnax
455  IF (abs(rat-1d0).GT.1d-8.AND.np.EQ.-1) CALL sarea (eps,rat)
456  if (abs(rat-1d0).GT.1d-8.AND.np.GE.0) CALL surfch(np,eps,rat)
457  IF (abs(rat-1d0).GT.1d-8.AND.np.EQ.-2) CALL sareac (eps,rat)
458 C PRINT 8000, RAT
459  8000 FORMAT ('RAT=',f8.6)
460  IF(np.EQ.-1.AND.eps.GE.1d0) print 7000,eps
461  IF(np.EQ.-1.AND.eps.LT.1d0) print 7001,eps
462  IF(np.GE.0) print 7100,np,eps
463  IF(np.EQ.-2.AND.eps.GE.1d0) print 7150,eps
464  IF(np.EQ.-2.AND.eps.LT.1d0) print 7151,eps
465  print 7400, lam,mrr,mri
466  print 7200, ddelt
467  7000 FORMAT('RANDOMLY ORIENTED OBLATE SPHEROIDS, A/B=',f11.7)
468  7001 FORMAT('RANDOMLY ORIENTED PROLATE SPHEROIDS, A/B=',f11.7)
469  7100 FORMAT('RANDOMLY ORIENTED CHEBYSHEV PARTICLES, T',
470  & i1,'(',f5.2,')')
471  7150 FORMAT('RANDOMLY ORIENTED OBLATE CYLINDERS, D/L=',f11.7)
472  7151 FORMAT('RANDOMLY ORIENTED PROLATE CYLINDERS, D/L=',f11.7)
473  7200 FORMAT ('ACCURACY OF COMPUTATIONS DDELT = ',d8.2)
474  7400 FORMAT('LAM=',f10.6,3x,'MRR=',d10.4,3x,'MRI=',d10.4)
475  ddelt=0.1d0*ddelt
476  DO 600 iax=1,npnax
477  axi=axmax-dax*dfloat(iax-1)
478  sz(iax) = axi
479  r1=0.89031d0*axi
480  r2=1.56538d0*axi
481  nk=int(axi*nkmax/axmax+2)
482  IF (nk.GT.1000) print 8001,nk
483  IF (nk.GT.1000) RETURN
484  IF (ndistr.EQ.3) CALL power (axi,b,r1,r2)
485  8001 FORMAT ('NK=',i4,' I.E., IS GREATER THAN 1000. ',
486  & 'EXECUTION TERMINATED.')
487  CALL gauss (nk,0,0,xg,wg)
488  z1=(r2-r1)*0.5d0
489  z2=(r1+r2)*0.5d0
490  z3=r1*0.5d0
491  IF (ndistr.EQ.5) GO TO 3
492  DO i=1,nk
493  xg1(i)=z1*xg(i)+z2
494  wg1(i)=wg(i)*z1
495  ENDDO
496  GO TO 4
497  3 DO i=1,nk
498  xg1(i)=z3*xg(i)+z3
499  wg1(i)=wg(i)*z3
500  ENDDO
501  DO i=nk+1,2*nk
502  ii=i-nk
503  xg1(i)=z1*xg(ii)+z2
504  wg1(i)=wg(ii)*z1
505  ENDDO
506  nk=nk*2
507  4 CALL distrb (nk,xg1,wg1,ndistr,axi,b,gam,r1,r2,
508  & reff(iax),veff(iax),p)
509  print 8002,r1,r2
510  8002 FORMAT('R1=',f10.6,' R2=',f10.6)
511  IF (abs(rat-1d0).LE.1d-6) print 8003, reff(iax),veff(iax)
512  IF (abs(rat-1d0).GT.1d-6) print 8004, reff(iax),veff(iax)
513  8003 FORMAT('EQUAL-VOLUME-SPHERE REFF=',f8.4,' VEFF=',f7.4)
514  8004 FORMAT('EQUAL-SURFACE-AREA-SPHERE REFF=',f8.4,
515  & ' VEFF=',f7.4)
516  print 7250,nk
517  7250 FORMAT('NUMBER OF GAUSSIAN QUADRATURE POINTS ',
518  & 'IN SIZE AVERAGING =',i4)
519  DO i=1,npl
520  alph1(i)=0d0
521  alph2(i)=0d0
522  alph3(i)=0d0
523  alph4(i)=0d0
524  bet1(i)=0d0
525  bet2(i)=0d0
526  alpha(1,i,iax)=0d0
527  alpha(2,i,iax)=0d0
528  alpha(3,i,iax)=0d0
529  alpha(4,i,iax)=0d0
530  beta(1,i,iax)=0d0
531  beta(2,i,iax)=0d0
532  ENDDO
533  cscat(iax)=0d0
534  cextin(iax)=0d0
535  l1max=0
536  DO 500 ink=1,nk
537  i=nk-ink+1
538  a=rat*xg1(i)
539  xev=2d0*p*a/lam
540  ixxx=xev+4.05d0*xev**0.333333d0
541  inm1=max0(4,ixxx)
542  IF (inm1.GE.npn1) THEN
543  print 7333, inm1, npn1, a
544  reff(iax)=-999.9
545  veff(iax)=-999.9
546  cextin(iax)=-999.9
547  cscat(iax)=-999.9
548  walb(iax)=-999.9
549  asymm(iax)=-999.9
550  fm(:,:,iax)=-999.9
551  GOTO 600
552  ENDIF
553  7333 FORMAT('INM1 (',i3,') is greater than NPN1 (',i3,') for radius ',d8.2)
554  qext1=0d0
555  qsca1=0d0
556  DO 50 nma=inm1,npn1
557  nmax=nma
558  mmax=1
559  ngauss=nmax*ndgs
560  IF (ngauss.GT.npng1) THEN
561  print 7340, ngauss, npng1, a
562  reff(iax)=-999.9
563  veff(iax)=-999.9
564  cextin(iax)=-999.9
565  cscat(iax)=-999.9
566  walb(iax)=-999.9
567  asymm(iax)=-999.9
568  fm(:,:,iax)=-999.9
569  GOTO 600
570  ENDIF
571  7340 FORMAT('NGAUSS (',i3,') is greater than NPNG1 (',i3,') for radius ',d8.2)
572  7334 FORMAT(' NMAX =', i3,' DC2 =',d8.2,' DC1 =',d8.2)
573  7335 FORMAT(' NMAX1 =', i3,' DC2 =',d8.2,' DC1 =',d8.2)
574  CALL const(ngauss,nmax,mmax,p,x,w,an,ann,s,ss,np,eps)
575  CALL vary(lam,mrr,mri,a,eps,np,ngauss,x,p,ppi,pir,pii,r,
576  & dr,ddr,drr,dri,nmax)
577  CALL tmatr0 (ngauss,x,w,an,ann,s,ss,ppi,pir,pii,r,dr,
578  & ddr,drr,dri,nmax,ncheck)
579  qext=0d0
580  qsca=0d0
581  DO n=1,nmax
582  n1=n+nmax
583  tr1nn=tr1(n,n)
584  ti1nn=ti1(n,n)
585  tr1nn1=tr1(n1,n1)
586  ti1nn1=ti1(n1,n1)
587  dn1=dfloat(2*n+1)
588  qsca=qsca+dn1*(tr1nn*tr1nn+ti1nn*ti1nn
589  & +tr1nn1*tr1nn1+ti1nn1*ti1nn1)
590  qext=qext+(tr1nn+tr1nn1)*dn1
591  ENDDO
592  dsca=abs((qsca1-qsca)/qsca)
593  dext=abs((qext1-qext)/qext)
594 C PRINT 7334, NMAX,DSCA,DEXT
595  qext1=qext
596  qsca1=qsca
597  nmin=dfloat(nmax)/2d0+1d0
598  DO 10 n=nmin,nmax
599  n1=n+nmax
600  tr1nn=tr1(n,n)
601  ti1nn=ti1(n,n)
602  tr1nn1=tr1(n1,n1)
603  ti1nn1=ti1(n1,n1)
604  dn1=dfloat(2*n+1)
605  dqsca=dn1*(tr1nn*tr1nn+ti1nn*ti1nn
606  & +tr1nn1*tr1nn1+ti1nn1*ti1nn1)
607  dqext=(tr1nn+tr1nn1)*dn1
608  dqsca=abs(dqsca/qsca)
609  dqext=abs(dqext/qext)
610  nmax1=n
611  IF (dqsca.LE.ddelt.AND.dqext.LE.ddelt) GO TO 12
612  10 CONTINUE
613  12 CONTINUE
614 C PRINT 7335, NMAX1,DQSCA,DQEXT
615  IF(dsca.LE.ddelt.AND.dext.LE.ddelt) GO TO 55
616  IF (nma.EQ.npn1) THEN
617  print 7333, nma, npn1, a
618  reff(iax)=-999.9
619  veff(iax)=-999.9
620  cextin(iax)=-999.9
621  cscat(iax)=-999.9
622  walb(iax)=-999.9
623  asymm(iax)=-999.9
624  fm(:,:,iax)=-999.9
625  GOTO 600
626  ENDIF
627  50 CONTINUE
628  55 nnnggg=ngauss+1
629  IF (ngauss.GE.npng1) THEN
630  print 7336, ngauss, npng1, a
631  print 7337, ng, dsca, dext
632  reff(iax)=-999.9
633  veff(iax)=-999.9
634  cextin(iax)=-999.9
635  cscat(iax)=-999.9
636  walb(iax)=-999.9
637  asymm(iax)=-999.9
638  fm(:,:,iax)=-999.9
639  GOTO 600
640  ENDIF
641  mmax=nmax1
642  DO 150 ngaus=nnnggg,npng1
643  ngauss=ngaus
644  nggg=2*ngauss
645  7336 FORMAT('NGAUSS (',i3,') is greater than NPNG1 (',i3,') for size ',d8.2)
646  7337 FORMAT(' NG = ',i3,' DSCA = ',d8.2,' DEXT = ',d8.2)
647  CALL const(ngauss,nmax,mmax,p,x,w,an,ann,s,ss,np,eps)
648  CALL vary(lam,mrr,mri,a,eps,np,ngauss,x,p,ppi,pir,pii,r,
649  & dr,ddr,drr,dri,nmax)
650  CALL tmatr0 (ngauss,x,w,an,ann,s,ss,ppi,pir,pii,r,dr,
651  & ddr,drr,dri,nmax,ncheck)
652  qext=0d0
653  qsca=0d0
654  DO 104 n=1,nmax
655  n1=n+nmax
656  tr1nn=tr1(n,n)
657  ti1nn=ti1(n,n)
658  tr1nn1=tr1(n1,n1)
659  ti1nn1=ti1(n1,n1)
660  dn1=dfloat(2*n+1)
661  qsca=qsca+dn1*(tr1nn*tr1nn+ti1nn*ti1nn
662  & +tr1nn1*tr1nn1+ti1nn1*ti1nn1)
663  qext=qext+(tr1nn+tr1nn1)*dn1
664  104 CONTINUE
665  dsca=abs((qsca1-qsca)/qsca)
666  dext=abs((qext1-qext)/qext)
667 C PRINT 7337, NGGG,DSCA,DEXT
668  qext1=qext
669  qsca1=qsca
670  IF(dsca.LE.ddelt.AND.dext.LE.ddelt) GOTO 155
671  IF (ngaus.GE.npng1) THEN
672  print 7336, ngaus, npng1, a
673  print 7337, nggg, dsca, dext
674  reff(iax)=-999.9
675  veff(iax)=-999.9
676  cextin(iax)=-999.9
677  cscat(iax)=-999.9
678  walb(iax)=-999.9
679  asymm(iax)=-999.9
680  fm(:,:,iax)=-999.9
681  GOTO 600
682  ENDIF
683  150 CONTINUE
684  155 CONTINUE
685  qsca=0d0
686  qext=0d0
687  nnm=nmax*2
688  DO 204 n=1,nnm
689  qext=qext+tr1(n,n)
690  204 CONTINUE
691  IF (nmax1.GE.npn4) THEN
692  print 7550, nmax1, npn4, a
693  reff(iax)=-999.9
694  veff(iax)=-999.9
695  cextin(iax)=-999.9
696  cscat(iax)=-999.9
697  walb(iax)=-999.9
698  asymm(iax)=-999.9
699  fm(:,:,iax)=-999.9
700  GOTO 600
701  ENDIF
702  7550 FORMAT ('NMAX1 (',i3,') is greater than NPN4 (',i3,') for radius ',d8.2)
703  IF (nmax1.GT.npn4) RETURN
704  DO 213 n2=1,nmax1
705  nn2=n2+nmax
706  DO 213 n1=1,nmax1
707  nn1=n1+nmax
708  zz1=tr1(n1,n2)
709  rt11(1,n1,n2)=zz1
710  zz2=ti1(n1,n2)
711  it11(1,n1,n2)=zz2
712  zz3=tr1(n1,nn2)
713  rt12(1,n1,n2)=zz3
714  zz4=ti1(n1,nn2)
715  it12(1,n1,n2)=zz4
716  zz5=tr1(nn1,n2)
717  rt21(1,n1,n2)=zz5
718  zz6=ti1(nn1,n2)
719  it21(1,n1,n2)=zz6
720  zz7=tr1(nn1,nn2)
721  rt22(1,n1,n2)=zz7
722  zz8=ti1(nn1,nn2)
723  it22(1,n1,n2)=zz8
724  qsca=qsca+zz1*zz1+zz2*zz2+zz3*zz3+zz4*zz4
725  & +zz5*zz5+zz6*zz6+zz7*zz7+zz8*zz8
726  213 CONTINUE
727 C PRINT 7800,0,ABS(QEXT),QSCA,NMAX
728  DO 220 m=1,nmax1
729  CALL tmatr(m,ngauss,x,w,an,ann,s,ss,ppi,pir,pii,r,dr,
730  & ddr,drr,dri,nmax,ncheck)
731  nm=nmax-m+1
732  nm1=nmax1-m+1
733  m1=m+1
734  qsc=0d0
735  DO 214 n2=1,nm1
736  nn2=n2+m-1
737  n22=n2+nm
738  DO 214 n1=1,nm1
739  nn1=n1+m-1
740  n11=n1+nm
741  zz1=tr1(n1,n2)
742  rt11(m1,nn1,nn2)=zz1
743  zz2=ti1(n1,n2)
744  it11(m1,nn1,nn2)=zz2
745  zz3=tr1(n1,n22)
746  rt12(m1,nn1,nn2)=zz3
747  zz4=ti1(n1,n22)
748  it12(m1,nn1,nn2)=zz4
749  zz5=tr1(n11,n2)
750  rt21(m1,nn1,nn2)=zz5
751  zz6=ti1(n11,n2)
752  it21(m1,nn1,nn2)=zz6
753  zz7=tr1(n11,n22)
754  rt22(m1,nn1,nn2)=zz7
755  zz8=ti1(n11,n22)
756  it22(m1,nn1,nn2)=zz8
757  qsc=qsc+(zz1*zz1+zz2*zz2+zz3*zz3+zz4*zz4
758  & +zz5*zz5+zz6*zz6+zz7*zz7+zz8*zz8)*2d0
759  214 CONTINUE
760  nnm=2*nm
761  qxt=0d0
762  DO 215 n=1,nnm
763  qxt=qxt+tr1(n,n)*2d0
764  215 CONTINUE
765  qsca=qsca+qsc
766  qext=qext+qxt
767 C PRINT 7800,M,ABS(QXT),QSC,NMAX
768  7800 FORMAT(' m = ',i3,' qxt = ',d12.6,' qsc = ',d12.6,
769  & ' nmax =',i3)
770  220 CONTINUE
771  coeff1=lam*lam*0.5d0/p
772  csca=qsca*coeff1
773  cext=-qext*coeff1
774 C PRINT 7880, NMAX,NMAX1
775  7880 FORMAT ('nmax=',i3,' nmax1=',i3)
776  CALL gsp (nmax1,csca,lam,al1,al2,al3,al4,be1,be2,lmax)
777  l1m=lmax+1
778  l1max=max(l1max,l1m)
779  wgii=wg1(i)
780  wgi=wgii*csca
781  DO 250 l1=1,l1m
782  alph1(l1)=alph1(l1)+al1(l1)*wgi
783  alph2(l1)=alph2(l1)+al2(l1)*wgi
784  alph3(l1)=alph3(l1)+al3(l1)*wgi
785  alph4(l1)=alph4(l1)+al4(l1)*wgi
786  bet1(l1)=bet1(l1)+be1(l1)*wgi
787  bet2(l1)=bet2(l1)+be2(l1)*wgi
788  250 CONTINUE
789  cscat(iax)=cscat(iax)+wgi
790  cextin(iax)=cextin(iax)+cext*wgii
791 C PRINT 6070, I,NMAX,NMAX1,NGAUSS
792  6070 FORMAT(4i6)
793  500 CONTINUE
794  DO 510 l1=1,l1max
795  alph1(l1)=alph1(l1)/cscat(iax)
796  alph2(l1)=alph2(l1)/cscat(iax)
797  alph3(l1)=alph3(l1)/cscat(iax)
798  alph4(l1)=alph4(l1)/cscat(iax)
799  bet1(l1)=bet1(l1)/cscat(iax)
800  bet2(l1)=bet2(l1)/cscat(iax)
801  alpha(1,l1,iax)=alph1(l1)
802  alpha(2,l1,iax)=alph2(l1)
803  alpha(3,l1,iax)=alph3(l1)
804  alpha(4,l1,iax)=alph4(l1)
805  beta(1,l1,iax)=bet1(l1)
806  beta(2,l1,iax)=bet2(l1)
807  510 CONTINUE
808  walb(iax)=cscat(iax)/cextin(iax)
809  CALL hovenr(l1max,alph1,alph2,alph3,alph4,bet1,bet2)
810  asymm(iax)=alph1(2)/3d0
811  print 9100,cextin(iax),cscat(iax),walb(iax),asymm(iax)
812  9100 FORMAT('CEXT=',d12.6,2x,'CSCA=',d12.6,2x,
813  & 2x,'W=',d12.6,2x,'<COS>=',d12.6)
814  IF (walb(iax).GT.1d0) print 9111
815  9111 FORMAT ('WARNING: W IS GREATER THAN 1')
816  WRITE (10,580) walb(iax),l1max
817  DO l=1,l1max
818  WRITE (10,575) alph1(l),alph2(l),alph3(l),alph4(l),
819  & bet1(l),bet2(l)
820  ENDDO
821  575 FORMAT(6d14.7)
822  580 FORMAT(d14.8,i8)
823  lmax=l1max-1
824  CALL matr (alph1,alph2,alph3,alph4,bet1,bet2,lmax,npna,npnax,iax,fm)
825  600 CONTINUE
826  itime=mclock()
827  time=dfloat(itime)/6000d0
828  print 1001,time
829  1001 FORMAT (' time =',f8.2,' min')
830  RETURN
831  END
832 
833 C**********************************************************************
834 C *
835 C INPUT PARAMETERS: *
836 C *
837 C NG = 2*NGAUSS - number of quadrature points on the *
838 C interval (-1,1). NGAUSS.LE.NPNG1 *
839 C NMAX,MMAX - maximum dimensions of the arrays. NMAX.LE.NPN1 *
840 C MMAX.LE.NPN1 *
841 C P - pi *
842 C *
843 C OUTPUT PARAMETERS: *
844 C *
845 C X,W - points and weights of the quadrature formula *
846 C AN(N) = n*(n+1) *
847 C ANN(N1,N2) = (1/2)*sqrt((2*n1+1)*(2*n2+1)/(n1*(n1+1)*n2*(n2+1))) *
848 C S(I)=1/sin(arccos(x(i))) *
849 C SS(I)=S(I)**2 *
850 C *
851 C**********************************************************************
852 
853  SUBROUTINE const (NGAUSS,NMAX,MMAX,P,X,W,AN,ANN,S,SS,NP,EPS)
854  IMPLICIT REAL*8 (a-h,o-z)
855  include 'tm_tmd.par.f'
856  real*8 x(npng2),w(npng2),x1(npng1),w1(npng1),
857  * x2(npng1),w2(npng1),
858  * s(npng2),ss(npng2),
859  * an(npn1),ann(npn1,npn1),dd(npn1)
860 
861  DO 10 n=1,nmax
862  nn=n*(n+1)
863  an(n)=dfloat(nn)
864  d=dsqrt(dfloat(2*n+1)/dfloat(nn))
865  dd(n)=d
866  DO 10 n1=1,n
867  ddd=d*dd(n1)*0.5d0
868  ann(n,n1)=ddd
869  ann(n1,n)=ddd
870  10 CONTINUE
871  ng=2*ngauss
872  IF (np.EQ.-2) GO TO 11
873  CALL gauss(ng,0,0,x,w)
874  GO TO 19
875  11 ng1=dfloat(ngauss)/2d0
876  ng2=ngauss-ng1
877  xx=-dcos(datan(eps))
878  CALL gauss(ng1,0,0,x1,w1)
879  CALL gauss(ng2,0,0,x2,w2)
880  DO 12 i=1,ng1
881  w(i)=0.5d0*(xx+1d0)*w1(i)
882  x(i)=0.5d0*(xx+1d0)*x1(i)+0.5d0*(xx-1d0)
883  12 CONTINUE
884  DO 14 i=1,ng2
885  w(i+ng1)=-0.5d0*xx*w2(i)
886  x(i+ng1)=-0.5d0*xx*x2(i)+0.5d0*xx
887  14 CONTINUE
888  DO 16 i=1,ngauss
889  w(ng-i+1)=w(i)
890  x(ng-i+1)=-x(i)
891  16 CONTINUE
892  19 DO 20 i=1,ngauss
893  y=x(i)
894  y=1d0/(1d0-y*y)
895  ss(i)=y
896  ss(ng-i+1)=y
897  y=dsqrt(y)
898  s(i)=y
899  s(ng-i+1)=y
900  20 CONTINUE
901  RETURN
902  END
903 
904 C**********************************************************************
905 C *
906 C INPUT PARAMETERS: *
907 C *
908 C LAM - wavelength of light *
909 C MRR and MRI - real and imaginary parts of the refractive index *
910 C A,EPS,NP - specify shape of the particle *
911 C (see subroutines RSP1, RSP2, and RSP3) *
912 C NG = NGAUSS*2 - number of gaussian quadrature points on the *
913 C interval (-1,1) *
914 C X - gaussian division points *
915 C P - pi *
916 C *
917 C OUTPUT INFORMATION: *
918 C *
919 C PPI = PI**2 , where PI = (2*P)/LAM (wavenumber) *
920 C PIR = PPI*MRR *
921 C PII = PPI*MRI *
922 C R and DR - see subroutines RSP1, RSP2, and RSP3 *
923 C DDR=1/(PI*SQRT(R)) *
924 C DRR+I*DRI=DDR/(MRR+I*MRI) *
925 C NMAX - dimension of T(m)-matrices *
926 C arrays J,Y,JR,JI,DJ,DY,DJR,DJI are transferred through *
927 C COMMON /CBESS/ - see subroutine BESS *
928 C *
929 C**********************************************************************
930 
931  SUBROUTINE vary (LAM,MRR,MRI,A,EPS,NP,NGAUSS,X,P,PPI,PIR,PII,
932  * R,DR,DDR,DRR,DRI,NMAX)
933  include 'tm_tmd.par.f'
934  IMPLICIT REAL*8 (a-h,o-z)
935  real*8 x(npng2),r(npng2),dr(npng2),mrr,mri,lam,
936  * z(npng2),zr(npng2),zi(npng2),
937  * j(npng2,npn1),y(npng2,npn1),jr(npng2,npn1),
938  * ji(npng2,npn1),dj(npng2,npn1),
939  * djr(npng2,npn1),dji(npng2,npn1),ddr(npng2),
940  * drr(npng2),dri(npng2),
941  * dy(npng2,npn1)
942  COMMON /cbess/ j,y,jr,ji,dj,dy,djr,dji
943  ng=ngauss*2
944  IF (np.EQ.-1) CALL rsp1(x,ng,ngauss,a,eps,np,r,dr)
945  IF (np.GE.0) CALL rsp2(x,ng,a,eps,np,r,dr)
946  IF (np.EQ.-2) CALL rsp3(x,ng,ngauss,a,eps,r,dr)
947  pi=p*2d0/lam
948  ppi=pi*pi
949  pir=ppi*mrr
950  pii=ppi*mri
951  v=1d0/(mrr*mrr+mri*mri)
952  prr=mrr*v
953  pri=-mri*v
954  ta=0d0
955  DO 10 i=1,ng
956  vv=dsqrt(r(i))
957  v=vv*pi
958  ta=max(ta,v)
959  vv=1d0/v
960  ddr(i)=vv
961  drr(i)=prr*vv
962  dri(i)=pri*vv
963  v1=v*mrr
964  v2=v*mri
965  z(i)=v
966  zr(i)=v1
967  zi(i)=v2
968  10 CONTINUE
969  IF (nmax.GT.npn1) print 9000,nmax,npn1
970  IF (nmax.GT.npn1) stop
971  9000 FORMAT(' NMAX = ',i2,', i.e., greater than ',i3)
972  tb=ta*dsqrt(mrr*mrr+mri*mri)
973  tb=dmax1(tb,dfloat(nmax))
974  nnmax1=1.2d0*dsqrt(dmax1(ta,dfloat(nmax)))+3d0
975  nnmax2=(tb+4d0*(tb**0.33333d0)+1.2d0*dsqrt(tb))
976  nnmax2=nnmax2-nmax+5
977  CALL bess(z,zr,zi,ng,nmax,nnmax1,nnmax2)
978  RETURN
979  END
980 
981 C**********************************************************************
982 C *
983 C Calculation of the functions R(I)=r(y(I))**2 and *
984 C DR(I)=((d/dy)r(y))/r(y) and horizontal semi-axis A *
985 C for a spheroid specified by the parameters REV (equal-volume- *
986 C sphere radius) and EPS=A/B (ratio of the semi-axes). *
987 C Y(I)=arccos(X(I)) *
988 C 1.LE.I.LE.NG *
989 C X - arguments *
990 C *
991 C**********************************************************************
992 
993  SUBROUTINE rsp1 (X,NG,NGAUSS,REV,EPS,NP,R,DR)
994  IMPLICIT REAL*8 (a-h,o-z)
995  real*8 x(ng),r(ng),dr(ng)
996  a=rev*eps**(1d0/3d0)
997  aa=a*a
998  ee=eps*eps
999  ee1=ee-1d0
1000  DO 50 i=1,ngauss
1001  c=x(i)
1002  cc=c*c
1003  ss=1d0-cc
1004  s=dsqrt(ss)
1005  rr=1d0/(ss+ee*cc)
1006  r(i)=aa*rr
1007  r(ng-i+1)=r(i)
1008  dr(i)=rr*c*s*ee1
1009  dr(ng-i+1)=-dr(i)
1010  50 CONTINUE
1011  RETURN
1012  END
1013 
1014 C**********************************************************************
1015 C *
1016 C Calculation of the functions R(I)=r(y(i))**2 and *
1017 C DR(I)=((d/dy)r(y))/r(y) and parameter R0 for a Chebyshev *
1018 C particle specified by the parameters REV (equal-volume-sphere *
1019 C radius), EPS, and N. *
1020 C Y(I)=arccos(X(I)) *
1021 C 1.LE.I.LE.NG *
1022 C X - arguments *
1023 C *
1024 C**********************************************************************
1025 
1026  SUBROUTINE rsp2 (X,NG,REV,EPS,N,R,DR)
1027  IMPLICIT REAL*8 (a-h,o-z)
1028  real*8 x(ng),r(ng),dr(ng)
1029  dnp=dfloat(n)
1030  dn=dnp*dnp
1031  dn4=dn*4d0
1032  ep=eps*eps
1033  a=1d0+1.5d0*ep*(dn4-2d0)/(dn4-1d0)
1034  i=(dnp+0.1d0)*0.5d0
1035  i=2*i
1036  IF (i.EQ.n) a=a-3d0*eps*(1d0+0.25d0*ep)/
1037  * (dn-1d0)-0.25d0*ep*eps/(9d0*dn-1d0)
1038  r0=rev*a**(-1d0/3d0)
1039  DO 50 i=1,ng
1040  xi=dacos(x(i))*dnp
1041  ri=r0*(1d0+eps*dcos(xi))
1042  r(i)=ri*ri
1043  dr(i)=-r0*eps*dnp*dsin(xi)/ri
1044  50 CONTINUE
1045  RETURN
1046  END
1047 
1048 C**********************************************************************
1049 C *
1050 C Calculation of the functions R(I)=R(Y(I))**2 and *
1051 C DR(I)=((d/dy)r(y))/r(y) *
1052 C for a cylinder specified by the parameters REV (equal-volume- *
1053 C sphere radius) and EPS=A/H (ratio of radius to semi-length) *
1054 C Y(I)=arccos(X(I)) *
1055 C 1.LE.I.LE.NG *
1056 C X - arguments *
1057 C *
1058 C**********************************************************************
1059 
1060  SUBROUTINE rsp3 (X,NG,NGAUSS,REV,EPS,R,DR)
1061  IMPLICIT REAL*8 (a-h,o-z)
1062  real*8 x(ng),r(ng),dr(ng)
1063  h=rev*( (2d0/(3d0*eps*eps))**(1d0/3d0) )
1064  a=h*eps
1065  DO 50 i=1,ngauss
1066  co=-x(i)
1067  si=dsqrt(1d0-co*co)
1068  IF (si/co.GT.a/h) GO TO 20
1069  rad=h/co
1070  rthet=h*si/(co*co)
1071  GO TO 30
1072  20 rad=a/si
1073  rthet=-a*co/(si*si)
1074  30 r(i)=rad*rad
1075  r(ng-i+1)=r(i)
1076  dr(i)=-rthet/rad
1077  dr(ng-i+1)=-dr(i)
1078  50 CONTINUE
1079  RETURN
1080  END
1081 
1082 C**********************************************************************
1083 C *
1084 C Calculation of spherical Bessel functions of the first kind *
1085 C J(I,N) = j_n(x) and second kind Y(I,N) = y_n(x) *
1086 C of real-valued argument X(I) and first kind JR(I,N)+I*JI(I,N) = *
1087 C = j_n(z) of complex argument Z(I)=XR(I)+I*XI(I), as well as *
1088 C the functions *
1089 C *
1090 C DJ(I,N) = (1/x)(d/dx)(x*j_n(x)) , *
1091 C DY(I,N) = (1/x)(d/dx)(x*y_n(x)) , *
1092 C DJR(I,N) = Re ((1/z)(d/dz)(z*j_n(x)) , *
1093 C DJI(I,N) = Im ((1/z)(d/dz)(z*j_n(x)) . *
1094 C *
1095 C 1.LE.N.LE.NMAX *
1096 C NMAX.LE.NPN1 *
1097 C X,XR,XI - arguments *
1098 C 1.LE.I.LE.NG *
1099 C Arrays J,Y,JR,JI,DJ,DY,DJR,DJI are in *
1100 C COMMON /CBESS/ *
1101 C Parameters NNMAX1 and NMAX2 determine computational accuracy *
1102 C *
1103 C**********************************************************************
1104 
1105  SUBROUTINE bess (X,XR,XI,NG,NMAX,NNMAX1,NNMAX2)
1106  include 'tm_tmd.par.f'
1107  IMPLICIT REAL*8 (a-h,o-z)
1108  real*8 x(ng),xr(ng),xi(ng),
1109  * j(npng2,npn1),y(npng2,npn1),jr(npng2,npn1),
1110  * ji(npng2,npn1),dj(npng2,npn1),dy(npng2,npn1),
1111  * djr(npng2,npn1),dji(npng2,npn1),
1112  * aj(npn1),ay(npn1),ajr(npn1),aji(npn1),
1113  * adj(npn1),ady(npn1),adjr(npn1),
1114  * adji(npn1)
1115  COMMON /cbess/ j,y,jr,ji,dj,dy,djr,dji
1116 
1117  DO 10 i=1,ng
1118  xx=x(i)
1119  CALL rjb(xx,aj,adj,nmax,nnmax1)
1120  CALL ryb(xx,ay,ady,nmax)
1121  yr=xr(i)
1122  yi=xi(i)
1123  CALL cjb(yr,yi,ajr,aji,adjr,adji,nmax,nnmax2)
1124  DO 10 n=1,nmax
1125  j(i,n)=aj(n)
1126  y(i,n)=ay(n)
1127  jr(i,n)=ajr(n)
1128  ji(i,n)=aji(n)
1129  dj(i,n)=adj(n)
1130  dy(i,n)=ady(n)
1131  djr(i,n)=adjr(n)
1132  dji(i,n)=adji(n)
1133  10 CONTINUE
1134  RETURN
1135  END
1136 
1137 C**********************************************************************
1138 C *
1139 C Calculation of spherical Bessel functions of the first kind j *
1140 C of real-valued argument x of orders from 1 to NMAX by using *
1141 C backward recursion. Parametr NNMAX determines numerical accuracy. *
1142 C U - function (1/x)(d/dx)(x*j(x)) *
1143 C *
1144 C**********************************************************************
1145 
1146  SUBROUTINE rjb(X,Y,U,NMAX,NNMAX)
1147  IMPLICIT REAL*8 (a-h,o-z)
1148  real*8 y(nmax),u(nmax),z(800)
1149  l=nmax+nnmax
1150  xx=1d0/x
1151  z(l)=1d0/(dfloat(2*l+1)*xx)
1152  l1=l-1
1153  DO 5 i=1,l1
1154  i1=l-i
1155  z(i1)=1d0/(dfloat(2*i1+1)*xx-z(i1+1))
1156  5 CONTINUE
1157  z0=1d0/(xx-z(1))
1158  y0=z0*dcos(x)*xx
1159  y1=y0*z(1)
1160  u(1)=y0-y1*xx
1161  y(1)=y1
1162  DO 10 i=2,nmax
1163  yi1=y(i-1)
1164  yi=yi1*z(i)
1165  u(i)=yi1-dfloat(i)*yi*xx
1166  y(i)=yi
1167  10 CONTINUE
1168  RETURN
1169  END
1170 
1171 C**********************************************************************
1172 C *
1173 C Calculation of spherical Bessel functions of the second kind y *
1174 C of real-valued argument x of orders from 1 to NMAX by using upward*
1175 C recursion. V - function (1/x)(d/dx)(x*y(x)) *
1176 C *
1177 C**********************************************************************
1178 
1179  SUBROUTINE ryb(X,Y,V,NMAX)
1180  IMPLICIT REAL*8 (a-h,o-z)
1181  real*8 y(nmax),v(nmax)
1182  c=dcos(x)
1183  s=dsin(x)
1184  x1=1d0/x
1185  x2=x1*x1
1186  x3=x2*x1
1187  y1=-c*x2-s*x1
1188  y(1)=y1
1189  y(2)=(-3d0*x3+x1)*c-3d0*x2*s
1190  nmax1=nmax-1
1191  DO 5 i=2,nmax1
1192  5 y(i+1)=dfloat(2*i+1)*x1*y(i)-y(i-1)
1193  v(1)=-x1*(c+y1)
1194  DO 10 i=2,nmax
1195  10 v(i)=y(i-1)-dfloat(i)*x1*y(i)
1196  RETURN
1197  END
1198 
1199 C**********************************************************************
1200 C *
1201 C Calculation of spherical Bessel functions of the first kind *
1202 C j=JR+I*JI of complex argument x=XR+I*XI of orders from 1 to NMAX *
1203 C by using backward recursion. Parametr NNMAX determines numerical *
1204 C accuracy. U=UR+I*UI - function (1/x)(d/dx)(x*j(x)) *
1205 C *
1206 C**********************************************************************
1207 
1208  SUBROUTINE cjb (XR,XI,YR,YI,UR,UI,NMAX,NNMAX)
1209  include 'tm_tmd.par.f'
1210  IMPLICIT REAL*8 (a-h,o-z)
1211  real*8 yr(nmax),yi(nmax),ur(nmax),ui(nmax)
1212  real*8 cyr(npn1),cyi(npn1),czr(1200),czi(1200),
1213  * cur(npn1),cui(npn1)
1214  l=nmax+nnmax
1215  xrxi=1d0/(xr*xr+xi*xi)
1216  cxxr=xr*xrxi
1217  cxxi=-xi*xrxi
1218  qf=1d0/dfloat(2*l+1)
1219  czr(l)=xr*qf
1220  czi(l)=xi*qf
1221  l1=l-1
1222  DO i=1,l1
1223  i1=l-i
1224  qf=dfloat(2*i1+1)
1225  ar=qf*cxxr-czr(i1+1)
1226  ai=qf*cxxi-czi(i1+1)
1227  ari=1d0/(ar*ar+ai*ai)
1228  czr(i1)=ar*ari
1229  czi(i1)=-ai*ari
1230  ENDDO
1231  ar=cxxr-czr(1)
1232  ai=cxxi-czi(1)
1233  ari=1d0/(ar*ar+ai*ai)
1234  cz0r=ar*ari
1235  cz0i=-ai*ari
1236  cr=dcos(xr)*dcosh(xi)
1237  ci=-dsin(xr)*dsinh(xi)
1238  ar=cz0r*cr-cz0i*ci
1239  ai=cz0i*cr+cz0r*ci
1240  cy0r=ar*cxxr-ai*cxxi
1241  cy0i=ai*cxxr+ar*cxxi
1242  cy1r=cy0r*czr(1)-cy0i*czi(1)
1243  cy1i=cy0i*czr(1)+cy0r*czi(1)
1244  ar=cy1r*cxxr-cy1i*cxxi
1245  ai=cy1i*cxxr+cy1r*cxxi
1246  cu1r=cy0r-ar
1247  cu1i=cy0i-ai
1248  cyr(1)=cy1r
1249  cyi(1)=cy1i
1250  cur(1)=cu1r
1251  cui(1)=cu1i
1252  yr(1)=cy1r
1253  yi(1)=cy1i
1254  ur(1)=cu1r
1255  ui(1)=cu1i
1256  DO i=2,nmax
1257  qi=dfloat(i)
1258  cyi1r=cyr(i-1)
1259  cyi1i=cyi(i-1)
1260  cyir=cyi1r*czr(i)-cyi1i*czi(i)
1261  cyii=cyi1i*czr(i)+cyi1r*czi(i)
1262  ar=cyir*cxxr-cyii*cxxi
1263  ai=cyii*cxxr+cyir*cxxi
1264  cuir=cyi1r-qi*ar
1265  cuii=cyi1i-qi*ai
1266  cyr(i)=cyir
1267  cyi(i)=cyii
1268  cur(i)=cuir
1269  cui(i)=cuii
1270  yr(i)=cyir
1271  yi(i)=cyii
1272  ur(i)=cuir
1273  ui(i)=cuii
1274  ENDDO
1275  RETURN
1276  END
1277 
1278 C**********************************************************************
1279 C *
1280 C calculation of the T(0) matrix for an axially symmetric particle *
1281 C *
1282 C Output information: *
1283 C *
1284 C Arrays TR1 and TI1 (real and imaginary parts of the *
1285 C T(0) matrix) are in COMMON /CT/ *
1286 C *
1287 C**********************************************************************
1288 
1289  SUBROUTINE tmatr0 (NGAUSS,X,W,AN,ANN,S,SS,PPI,PIR,PII,R,DR,DDR,
1290  * DRR,DRI,NMAX,NCHECK)
1291  include 'tm_tmd.par.f'
1292  IMPLICIT REAL*8 (a-h,o-z)
1293  real*8 x(npng2),w(npng2),an(npn1),s(npng2),ss(npng2),
1294  * r(npng2),dr(npng2),sig(npn2),
1295  * j(npng2,npn1),y(npng2,npn1),
1296  * jr(npng2,npn1),ji(npng2,npn1),dj(npng2,npn1),
1297  * dy(npng2,npn1),djr(npng2,npn1),
1298  * dji(npng2,npn1),ddr(npng2),drr(npng2),
1299  * d1(npng2,npn1),d2(npng2,npn1),
1300  * dri(npng2),ds(npng2),dss(npng2),rr(npng2),
1301  * dv1(npn1),dv2(npn1)
1302 
1303  real*8 r11(npn1,npn1),r12(npn1,npn1),
1304  * r21(npn1,npn1),r22(npn1,npn1),
1305  * i11(npn1,npn1),i12(npn1,npn1),
1306  * i21(npn1,npn1),i22(npn1,npn1),
1307  * rg11(npn1,npn1),rg12(npn1,npn1),
1308  * rg21(npn1,npn1),rg22(npn1,npn1),
1309  * ig11(npn1,npn1),ig12(npn1,npn1),
1310  * ig21(npn1,npn1),ig22(npn1,npn1),
1311  * ann(npn1,npn1),
1312  * qr(npn2,npn2),qi(npn2,npn2),
1313  * rgqr(npn2,npn2),rgqi(npn2,npn2),
1314  * tqr(npn2,npn2),tqi(npn2,npn2),
1315  * trgqr(npn2,npn2),trgqi(npn2,npn2)
1316  real*8 tr1(npn2,npn2),ti1(npn2,npn2)
1317  real*4 plus(npn6*npn4*npn4*8)
1318  COMMON /tmat/ plus,
1319  & r11,r12,r21,r22,i11,i12,i21,i22,rg11,rg12,rg21,rg22,
1320  & ig11,ig12,ig21,ig22
1321  COMMON /cbess/ j,y,jr,ji,dj,dy,djr,dji
1322  COMMON /ct/ tr1,ti1
1323  COMMON /ctt/ qr,qi,rgqr,rgqi
1324  mm1=1
1325  nnmax=nmax+nmax
1326  ng=2*ngauss
1327  ngss=ng
1328  factor=1d0
1329  IF (ncheck.EQ.1) THEN
1330  ngss=ngauss
1331  factor=2d0
1332  ELSE
1333  CONTINUE
1334  ENDIF
1335  si=1d0
1336  DO 5 n=1,nnmax
1337  si=-si
1338  sig(n)=si
1339  5 CONTINUE
1340  20 DO 25 i=1,ngauss
1341  i1=ngauss+i
1342  i2=ngauss-i+1
1343  CALL vig ( x(i1), nmax, 0, dv1, dv2)
1344  DO 25 n=1,nmax
1345  si=sig(n)
1346  dd1=dv1(n)
1347  dd2=dv2(n)
1348  d1(i1,n)=dd1
1349  d2(i1,n)=dd2
1350  d1(i2,n)=dd1*si
1351  d2(i2,n)=-dd2*si
1352  25 CONTINUE
1353  30 DO 40 i=1,ngss
1354  rr(i)=w(i)*r(i)
1355  40 CONTINUE
1356 
1357  DO 300 n1=mm1,nmax
1358  an1=an(n1)
1359  DO 300 n2=mm1,nmax
1360  an2=an(n2)
1361  ar12=0d0
1362  ar21=0d0
1363  ai12=0d0
1364  ai21=0d0
1365  gr12=0d0
1366  gr21=0d0
1367  gi12=0d0
1368  gi21=0d0
1369  IF (ncheck.EQ.1.AND.sig(n1+n2).LT.0d0) GO TO 205
1370 
1371  DO 200 i=1,ngss
1372  d1n1=d1(i,n1)
1373  d2n1=d2(i,n1)
1374  d1n2=d1(i,n2)
1375  d2n2=d2(i,n2)
1376  a12=d1n1*d2n2
1377  a21=d2n1*d1n2
1378  a22=d2n1*d2n2
1379  aa1=a12+a21
1380 
1381  qj1=j(i,n1)
1382  qy1=y(i,n1)
1383  qjr2=jr(i,n2)
1384  qji2=ji(i,n2)
1385  qdjr2=djr(i,n2)
1386  qdji2=dji(i,n2)
1387  qdj1=dj(i,n1)
1388  qdy1=dy(i,n1)
1389 
1390  c1r=qjr2*qj1
1391  c1i=qji2*qj1
1392  b1r=c1r-qji2*qy1
1393  b1i=c1i+qjr2*qy1
1394 
1395  c2r=qjr2*qdj1
1396  c2i=qji2*qdj1
1397  b2r=c2r-qji2*qdy1
1398  b2i=c2i+qjr2*qdy1
1399 
1400  ddri=ddr(i)
1401  c3r=ddri*c1r
1402  c3i=ddri*c1i
1403  b3r=ddri*b1r
1404  b3i=ddri*b1i
1405 
1406  c4r=qdjr2*qj1
1407  c4i=qdji2*qj1
1408  b4r=c4r-qdji2*qy1
1409  b4i=c4i+qdjr2*qy1
1410 
1411  drri=drr(i)
1412  drii=dri(i)
1413  c5r=c1r*drri-c1i*drii
1414  c5i=c1i*drri+c1r*drii
1415  b5r=b1r*drri-b1i*drii
1416  b5i=b1i*drri+b1r*drii
1417 
1418  uri=dr(i)
1419  rri=rr(i)
1420 
1421  f1=rri*a22
1422  f2=rri*uri*an1*a12
1423  ar12=ar12+f1*b2r+f2*b3r
1424  ai12=ai12+f1*b2i+f2*b3i
1425  gr12=gr12+f1*c2r+f2*c3r
1426  gi12=gi12+f1*c2i+f2*c3i
1427 
1428  f2=rri*uri*an2*a21
1429  ar21=ar21+f1*b4r+f2*b5r
1430  ai21=ai21+f1*b4i+f2*b5i
1431  gr21=gr21+f1*c4r+f2*c5r
1432  gi21=gi21+f1*c4i+f2*c5i
1433  200 CONTINUE
1434 
1435  205 an12=ann(n1,n2)*factor
1436  r12(n1,n2)=ar12*an12
1437  r21(n1,n2)=ar21*an12
1438  i12(n1,n2)=ai12*an12
1439  i21(n1,n2)=ai21*an12
1440  rg12(n1,n2)=gr12*an12
1441  rg21(n1,n2)=gr21*an12
1442  ig12(n1,n2)=gi12*an12
1443  ig21(n1,n2)=gi21*an12
1444  300 CONTINUE
1445 
1446  tpir=pir
1447  tpii=pii
1448  tppi=ppi
1449 
1450  nm=nmax
1451  DO 310 n1=mm1,nmax
1452  k1=n1-mm1+1
1453  kk1=k1+nm
1454  DO 310 n2=mm1,nmax
1455  k2=n2-mm1+1
1456  kk2=k2+nm
1457 
1458  tar12= i12(n1,n2)
1459  tai12=-r12(n1,n2)
1460  tgr12= ig12(n1,n2)
1461  tgi12=-rg12(n1,n2)
1462 
1463  tar21=-i21(n1,n2)
1464  tai21= r21(n1,n2)
1465  tgr21=-ig21(n1,n2)
1466  tgi21= rg21(n1,n2)
1467 
1468  tqr(k1,k2)=tpir*tar21-tpii*tai21+tppi*tar12
1469  tqi(k1,k2)=tpir*tai21+tpii*tar21+tppi*tai12
1470  trgqr(k1,k2)=tpir*tgr21-tpii*tgi21+tppi*tgr12
1471  trgqi(k1,k2)=tpir*tgi21+tpii*tgr21+tppi*tgi12
1472 
1473  tqr(k1,kk2)=0d0
1474  tqi(k1,kk2)=0d0
1475  trgqr(k1,kk2)=0d0
1476  trgqi(k1,kk2)=0d0
1477 
1478  tqr(kk1,k2)=0d0
1479  tqi(kk1,k2)=0d0
1480  trgqr(kk1,k2)=0d0
1481  trgqi(kk1,k2)=0d0
1482 
1483  tqr(kk1,kk2)=tpir*tar12-tpii*tai12+tppi*tar21
1484  tqi(kk1,kk2)=tpir*tai12+tpii*tar12+tppi*tai21
1485  trgqr(kk1,kk2)=tpir*tgr12-tpii*tgi12+tppi*tgr21
1486  trgqi(kk1,kk2)=tpir*tgi12+tpii*tgr12+tppi*tgi21
1487  310 CONTINUE
1488 
1489  nnmax=2*nm
1490  DO 320 n1=1,nnmax
1491  DO 320 n2=1,nnmax
1492  qr(n1,n2)=tqr(n1,n2)
1493  qi(n1,n2)=tqi(n1,n2)
1494  rgqr(n1,n2)=trgqr(n1,n2)
1495  rgqi(n1,n2)=trgqi(n1,n2)
1496  320 CONTINUE
1497  CALL tt(nmax,ncheck)
1498  RETURN
1499  END
1500 
1501 C**********************************************************************
1502 C *
1503 C Calculation of the T(M) matrix, M.GE.1, for an axially symmetric *
1504 C particle *
1505 C *
1506 C Input parameters: *
1507 C *
1508 C M.GE.1 *
1509 C NG = NGAUSS*2 - number of gaussian division points on the *
1510 C interval (-1,1) *
1511 C W - quadrature weights *
1512 C AN,ANN - see subroutine CONST *
1513 C S,SS - see subroutine CONST *
1514 C ARRAYS DV1,DV2,DV3,DV4 are in COMMON /DV/ - *
1515 C see subroutine DVIG *
1516 C PPI,PIR,PII - see subroutine VARY *
1517 C R J DR - see subroutines RSP1 and RSP2 *
1518 C DDR,DRR,DRI - see subroutine VARY *
1519 C NMAX - dimension of the T(M) matrix *
1520 C Arrays J,Y,JR,JI,DJ,DY,DJR,DJI are in *
1521 C COMMON /CBESS/ - see subroutine BESS *
1522 C *
1523 C Output parameters: *
1524 C *
1525 C Arrays TR1,TI1 (real and imaginary parts of the T(M) matrix) *
1526 C are in COMMON /CT/ *
1527 C *
1528 C**********************************************************************
1529 
1530  SUBROUTINE tmatr (M,NGAUSS,X,W,AN,ANN,S,SS,PPI,PIR,PII,R,DR,DDR,
1531  * DRR,DRI,NMAX,NCHECK)
1532  include 'tm_tmd.par.f'
1533  IMPLICIT REAL*8 (a-h,o-z)
1534  real*8 x(npng2),w(npng2),an(npn1),s(npng2),ss(npng2),
1535  * r(npng2),dr(npng2),sig(npn2),
1536  * j(npng2,npn1),y(npng2,npn1),
1537  * jr(npng2,npn1),ji(npng2,npn1),dj(npng2,npn1),
1538  * dy(npng2,npn1),djr(npng2,npn1),
1539  * dji(npng2,npn1),ddr(npng2),drr(npng2),
1540  * d1(npng2,npn1),d2(npng2,npn1),
1541  * dri(npng2),ds(npng2),dss(npng2),rr(npng2),
1542  * dv1(npn1),dv2(npn1)
1543 
1544  real*8 r11(npn1,npn1),r12(npn1,npn1),
1545  * r21(npn1,npn1),r22(npn1,npn1),
1546  * i11(npn1,npn1),i12(npn1,npn1),
1547  * i21(npn1,npn1),i22(npn1,npn1),
1548  * rg11(npn1,npn1),rg12(npn1,npn1),
1549  * rg21(npn1,npn1),rg22(npn1,npn1),
1550  * ig11(npn1,npn1),ig12(npn1,npn1),
1551  * ig21(npn1,npn1),ig22(npn1,npn1),
1552  * ann(npn1,npn1),
1553  * qr(npn2,npn2),qi(npn2,npn2),
1554  * rgqr(npn2,npn2),rgqi(npn2,npn2),
1555  * tqr(npn2,npn2),tqi(npn2,npn2),
1556  * trgqr(npn2,npn2),trgqi(npn2,npn2)
1557  real*8 tr1(npn2,npn2),ti1(npn2,npn2)
1558  real*4 plus(npn6*npn4*npn4*8)
1559  COMMON /tmat/ plus,
1560  & r11,r12,r21,r22,i11,i12,i21,i22,rg11,rg12,rg21,rg22,
1561  & ig11,ig12,ig21,ig22
1562  COMMON /cbess/ j,y,jr,ji,dj,dy,djr,dji
1563  COMMON /ct/ tr1,ti1
1564  COMMON /ctt/ qr,qi,rgqr,rgqi
1565  mm1=m
1566  qm=dfloat(m)
1567  qmm=qm*qm
1568  ng=2*ngauss
1569  ngss=ng
1570  factor=1d0
1571  IF (ncheck.EQ.1) THEN
1572  ngss=ngauss
1573  factor=2d0
1574  ELSE
1575  CONTINUE
1576  ENDIF
1577  si=1d0
1578  nm=nmax+nmax
1579  DO 5 n=1,nm
1580  si=-si
1581  sig(n)=si
1582  5 CONTINUE
1583  20 DO 25 i=1,ngauss
1584  i1=ngauss+i
1585  i2=ngauss-i+1
1586  CALL vig (x(i1),nmax,m,dv1,dv2)
1587  DO 25 n=1,nmax
1588  si=sig(n)
1589  dd1=dv1(n)
1590  dd2=dv2(n)
1591  d1(i1,n)=dd1
1592  d2(i1,n)=dd2
1593  d1(i2,n)=dd1*si
1594  d2(i2,n)=-dd2*si
1595  25 CONTINUE
1596  30 DO 40 i=1,ngss
1597  wr=w(i)*r(i)
1598  ds(i)=s(i)*qm*wr
1599  dss(i)=ss(i)*qmm
1600  rr(i)=wr
1601  40 CONTINUE
1602 
1603  DO 300 n1=mm1,nmax
1604  an1=an(n1)
1605  DO 300 n2=mm1,nmax
1606  an2=an(n2)
1607  ar11=0d0
1608  ar12=0d0
1609  ar21=0d0
1610  ar22=0d0
1611  ai11=0d0
1612  ai12=0d0
1613  ai21=0d0
1614  ai22=0d0
1615  gr11=0d0
1616  gr12=0d0
1617  gr21=0d0
1618  gr22=0d0
1619  gi11=0d0
1620  gi12=0d0
1621  gi21=0d0
1622  gi22=0d0
1623  si=sig(n1+n2)
1624 
1625  DO 200 i=1,ngss
1626  d1n1=d1(i,n1)
1627  d2n1=d2(i,n1)
1628  d1n2=d1(i,n2)
1629  d2n2=d2(i,n2)
1630  a11=d1n1*d1n2
1631  a12=d1n1*d2n2
1632  a21=d2n1*d1n2
1633  a22=d2n1*d2n2
1634  aa1=a12+a21
1635  aa2=a11*dss(i)+a22
1636  qj1=j(i,n1)
1637  qy1=y(i,n1)
1638  qjr2=jr(i,n2)
1639  qji2=ji(i,n2)
1640  qdjr2=djr(i,n2)
1641  qdji2=dji(i,n2)
1642  qdj1=dj(i,n1)
1643  qdy1=dy(i,n1)
1644 
1645  c1r=qjr2*qj1
1646  c1i=qji2*qj1
1647  b1r=c1r-qji2*qy1
1648  b1i=c1i+qjr2*qy1
1649 
1650  c2r=qjr2*qdj1
1651  c2i=qji2*qdj1
1652  b2r=c2r-qji2*qdy1
1653  b2i=c2i+qjr2*qdy1
1654 
1655  ddri=ddr(i)
1656  c3r=ddri*c1r
1657  c3i=ddri*c1i
1658  b3r=ddri*b1r
1659  b3i=ddri*b1i
1660 
1661  c4r=qdjr2*qj1
1662  c4i=qdji2*qj1
1663  b4r=c4r-qdji2*qy1
1664  b4i=c4i+qdjr2*qy1
1665 
1666  drri=drr(i)
1667  drii=dri(i)
1668  c5r=c1r*drri-c1i*drii
1669  c5i=c1i*drri+c1r*drii
1670  b5r=b1r*drri-b1i*drii
1671  b5i=b1i*drri+b1r*drii
1672 
1673  c6r=qdjr2*qdj1
1674  c6i=qdji2*qdj1
1675  b6r=c6r-qdji2*qdy1
1676  b6i=c6i+qdjr2*qdy1
1677 
1678  c7r=c4r*ddri
1679  c7i=c4i*ddri
1680  b7r=b4r*ddri
1681  b7i=b4i*ddri
1682 
1683  c8r=c2r*drri-c2i*drii
1684  c8i=c2i*drri+c2r*drii
1685  b8r=b2r*drri-b2i*drii
1686  b8i=b2i*drri+b2r*drii
1687 
1688  uri=dr(i)
1689  dsi=ds(i)
1690  dssi=dss(i)
1691  rri=rr(i)
1692 
1693  IF (ncheck.EQ.1.AND.si.GT.0d0) GO TO 150
1694 
1695  e1=dsi*aa1
1696  ar11=ar11+e1*b1r
1697  ai11=ai11+e1*b1i
1698  gr11=gr11+e1*c1r
1699  gi11=gi11+e1*c1i
1700  IF (ncheck.EQ.1) GO TO 160
1701 
1702  150 f1=rri*aa2
1703  f2=rri*uri*an1*a12
1704  ar12=ar12+f1*b2r+f2*b3r
1705  ai12=ai12+f1*b2i+f2*b3i
1706  gr12=gr12+f1*c2r+f2*c3r
1707  gi12=gi12+f1*c2i+f2*c3i
1708 
1709  f2=rri*uri*an2*a21
1710  ar21=ar21+f1*b4r+f2*b5r
1711  ai21=ai21+f1*b4i+f2*b5i
1712  gr21=gr21+f1*c4r+f2*c5r
1713  gi21=gi21+f1*c4i+f2*c5i
1714  IF (ncheck.EQ.1) GO TO 200
1715 
1716  160 e2=dsi*uri*a11
1717  e3=e2*an2
1718  e2=e2*an1
1719  ar22=ar22+e1*b6r+e2*b7r+e3*b8r
1720  ai22=ai22+e1*b6i+e2*b7i+e3*b8i
1721  gr22=gr22+e1*c6r+e2*c7r+e3*c8r
1722  gi22=gi22+e1*c6i+e2*c7i+e3*c8i
1723  200 CONTINUE
1724  an12=ann(n1,n2)*factor
1725  r11(n1,n2)=ar11*an12
1726  r12(n1,n2)=ar12*an12
1727  r21(n1,n2)=ar21*an12
1728  r22(n1,n2)=ar22*an12
1729  i11(n1,n2)=ai11*an12
1730  i12(n1,n2)=ai12*an12
1731  i21(n1,n2)=ai21*an12
1732  i22(n1,n2)=ai22*an12
1733  rg11(n1,n2)=gr11*an12
1734  rg12(n1,n2)=gr12*an12
1735  rg21(n1,n2)=gr21*an12
1736  rg22(n1,n2)=gr22*an12
1737  ig11(n1,n2)=gi11*an12
1738  ig12(n1,n2)=gi12*an12
1739  ig21(n1,n2)=gi21*an12
1740  ig22(n1,n2)=gi22*an12
1741 
1742  300 CONTINUE
1743  tpir=pir
1744  tpii=pii
1745  tppi=ppi
1746  nm=nmax-mm1+1
1747  DO 310 n1=mm1,nmax
1748  k1=n1-mm1+1
1749  kk1=k1+nm
1750  DO 310 n2=mm1,nmax
1751  k2=n2-mm1+1
1752  kk2=k2+nm
1753 
1754  tar11=-r11(n1,n2)
1755  tai11=-i11(n1,n2)
1756  tgr11=-rg11(n1,n2)
1757  tgi11=-ig11(n1,n2)
1758 
1759  tar12= i12(n1,n2)
1760  tai12=-r12(n1,n2)
1761  tgr12= ig12(n1,n2)
1762  tgi12=-rg12(n1,n2)
1763 
1764  tar21=-i21(n1,n2)
1765  tai21= r21(n1,n2)
1766  tgr21=-ig21(n1,n2)
1767  tgi21= rg21(n1,n2)
1768 
1769  tar22=-r22(n1,n2)
1770  tai22=-i22(n1,n2)
1771  tgr22=-rg22(n1,n2)
1772  tgi22=-ig22(n1,n2)
1773 
1774  tqr(k1,k2)=tpir*tar21-tpii*tai21+tppi*tar12
1775  tqi(k1,k2)=tpir*tai21+tpii*tar21+tppi*tai12
1776  trgqr(k1,k2)=tpir*tgr21-tpii*tgi21+tppi*tgr12
1777  trgqi(k1,k2)=tpir*tgi21+tpii*tgr21+tppi*tgi12
1778 
1779  tqr(k1,kk2)=tpir*tar11-tpii*tai11+tppi*tar22
1780  tqi(k1,kk2)=tpir*tai11+tpii*tar11+tppi*tai22
1781  trgqr(k1,kk2)=tpir*tgr11-tpii*tgi11+tppi*tgr22
1782  trgqi(k1,kk2)=tpir*tgi11+tpii*tgr11+tppi*tgi22
1783 
1784  tqr(kk1,k2)=tpir*tar22-tpii*tai22+tppi*tar11
1785  tqi(kk1,k2)=tpir*tai22+tpii*tar22+tppi*tai11
1786  trgqr(kk1,k2)=tpir*tgr22-tpii*tgi22+tppi*tgr11
1787  trgqi(kk1,k2)=tpir*tgi22+tpii*tgr22+tppi*tgi11
1788 
1789  tqr(kk1,kk2)=tpir*tar12-tpii*tai12+tppi*tar21
1790  tqi(kk1,kk2)=tpir*tai12+tpii*tar12+tppi*tai21
1791  trgqr(kk1,kk2)=tpir*tgr12-tpii*tgi12+tppi*tgr21
1792  trgqi(kk1,kk2)=tpir*tgi12+tpii*tgr12+tppi*tgi21
1793  310 CONTINUE
1794 
1795  nnmax=2*nm
1796  DO 320 n1=1,nnmax
1797  DO 320 n2=1,nnmax
1798  qr(n1,n2)=tqr(n1,n2)
1799  qi(n1,n2)=tqi(n1,n2)
1800  rgqr(n1,n2)=trgqr(n1,n2)
1801  rgqi(n1,n2)=trgqi(n1,n2)
1802  320 CONTINUE
1803 
1804  CALL tt(nm,ncheck)
1805 
1806  RETURN
1807  END
1808 
1809 C*****************************************************************
1810 c
1811 c Calculation of the functiONS
1812 c DV1(n)=dvig(0,m,n,arccos x)
1813 c and
1814 c DV2(n)=[d/d(arccos x)] dvig(0,m,n,arccos x)
1815 c 1.LE.N.LE.NMAX
1816 c 0.LE.x.LE.1
1817 
1818  SUBROUTINE vig (X, NMAX, M, DV1, DV2)
1819  include 'tm_tmd.par.f'
1820  IMPLICIT REAL*8 (a-h,o-z)
1821  real*8 dv1(npn1),dv2(npn1)
1822 
1823  a=1d0
1824  qs=dsqrt(1d0-x*x)
1825  qs1=1d0/qs
1826  DO n=1,nmax
1827  dv1(n)=0d0
1828  dv2(n)=0d0
1829  ENDDO
1830  IF (m.NE.0) GO TO 20
1831  d1=1d0
1832  d2=x
1833  DO n=1,nmax
1834  qn=dfloat(n)
1835  qn1=dfloat(n+1)
1836  qn2=dfloat(2*n+1)
1837  d3=(qn2*x*d2-qn*d1)/qn1
1838  der=qs1*(qn1*qn/qn2)*(-d1+d3)
1839  dv1(n)=d2
1840  dv2(n)=der
1841  d1=d2
1842  d2=d3
1843  ENDDO
1844  RETURN
1845  20 qmm=dfloat(m*m)
1846  DO i=1,m
1847  i2=i*2
1848  a=a*dsqrt(dfloat(i2-1)/dfloat(i2))*qs
1849  ENDDO
1850  d1=0d0
1851  d2=a
1852  DO n=m,nmax
1853  qn=dfloat(n)
1854  qn2=dfloat(2*n+1)
1855  qn1=dfloat(n+1)
1856  qnm=dsqrt(qn*qn-qmm)
1857  qnm1=dsqrt(qn1*qn1-qmm)
1858  d3=(qn2*x*d2-qnm*d1)/qnm1
1859  der=qs1*(-qn1*qnm*d1+qn*qnm1*d3)/qn2
1860  dv1(n)=d2
1861  dv2(n)=der
1862  d1=d2
1863  d2=d3
1864  ENDDO
1865  RETURN
1866  END
1867 
1868 C**********************************************************************
1869 C *
1870 C Calculation of the matrix T = - RG(Q) * (Q**(-1)) *
1871 C *
1872 C Input infortmation is in COMMON /CTT/ *
1873 C Output information is in COMMON /CT/ *
1874 C *
1875 C**********************************************************************
1876 
1877  SUBROUTINE tt(NMAX,NCHECK)
1878  include 'tm_tmd.par.f'
1879  IMPLICIT REAL*8 (a-h,o-z)
1880  real*8 f(npn2,npn2),b(npn2),work(npn2),
1881  * qr(npn2,npn2),qi(npn2,npn2),
1882  * rgqr(npn2,npn2),rgqi(npn2,npn2),
1883  * a(npn2,npn2),c(npn2,npn2),d(npn2,npn2),e(npn2,npn2)
1884  real*8 tr1(npn2,npn2),ti1(npn2,npn2)
1885  COMPLEX*16 ZQ(NPN2,NPN2),ZW(NPN2)
1886  INTEGER IPIV(NPN2),IPVT(NPN2)
1887  COMMON /CT/ TR1,TI1
1888  COMMON /CTT/ QR,QI,RGQR,RGQI
1889  ndim=npn2
1890  nnmax=2*nmax
1891 
1892 C Matrix inversion from LAPACK
1893 
1894  DO i=1,nnmax
1895  DO j=1,nnmax
1896  zq(i,j)=dcmplx(qr(i,j),qi(i,j))
1897  ENDDO
1898  ENDDO
1899  info=0
1900  CALL zgetrf(nnmax,nnmax,zq,npn2,ipiv,info)
1901  IF (info.NE.0) WRITE (6,1100) info
1902  CALL zgetri(nnmax,zq,npn2,ipiv,zw,npn2,info)
1903  IF (info.NE.0) WRITE (6,1100) info
1904 
1905  1100 FORMAT ('WARNING: info=', i2)
1906  DO i=1,nnmax
1907  DO j=1,nnmax
1908  tr=0d0
1909  ti=0d0
1910  DO k=1,nnmax
1911  arr=rgqr(i,k)
1912  ari=rgqi(i,k)
1913  ar=zq(k,j)
1914  ai=dimag(zq(k,j))
1915  tr=tr-arr*ar+ari*ai
1916  ti=ti-arr*ai-ari*ar
1917  ENDDO
1918  tr1(i,j)=tr
1919  ti1(i,j)=ti
1920  ENDDO
1921  ENDDO
1922  RETURN
1923  END
1924 
1925 C********************************************************************
1926 C *
1927 C Calculation of the expansion coefficients for the (I,Q,U,V) *
1928 C representation. *
1929 C *
1930 C Input parameters: *
1931 C *
1932 C LAM - wavelength of light *
1933 C CSCA - scattering cross section *
1934 C TR and TI - elements of the t-matrix. Transferred through *
1935 C COMMON /CTM/ *
1936 C NMAX - dimension of T(M) matrices *
1937 C *
1938 C Output infortmation: *
1939 C *
1940 C ALF1,...,ALF4,BET1,BET2 - expansion coefficients *
1941 C LMAX - number of coefficients minus 1 *
1942 C *
1943 C********************************************************************
1944 
1945  SUBROUTINE gsp(NMAX,CSCA,LAM,ALF1,ALF2,ALF3,ALF4,BET1,BET2,LMAX)
1946  include 'tm_tmd.par.f'
1947  IMPLICIT REAL*8 (a-b,d-h,o-z),COMPLEX*16 (C)
1948  REAL*8 LAM,SSIGN(900)
1949  REAL*8 CSCA,SSI(NPL),SSJ(NPN1),
1950  & alf1(npl),alf2(npl),alf3(npl),
1951  & alf4(npl),bet1(npl),bet2(npl),
1952  & tr1(npl1,npn4),tr2(npl1,npn4),
1953  & ti1(npl1,npn4),ti2(npl1,npn4),
1954  & g1(npl1,npn6),g2(npl1,npn6),
1955  & ar1(npn4),ar2(npn4),ai1(npn4),ai2(npn4),
1956  & fr(npn4,npn4),fi(npn4,npn4),ff(npn4,npn4)
1957  real*4 b1r(npl1,npl1,npn4),b1i(npl1,npl1,npn4),
1958  & b2r(npl1,npl1,npn4),b2i(npl1,npl1,npn4),
1959  & d1(npl1,npn4,npn4),d2(npl1,npn4,npn4),
1960  & d3(npl1,npn4,npn4),d4(npl1,npn4,npn4),
1961  & d5r(npl1,npn4,npn4),d5i(npl1,npn4,npn4),
1962  & plus1(npn6*npn4*npn4*8)
1963  real*4
1964  & tr11(npn6,npn4,npn4),tr12(npn6,npn4,npn4),
1965  & tr21(npn6,npn4,npn4),tr22(npn6,npn4,npn4),
1966  & ti11(npn6,npn4,npn4),ti12(npn6,npn4,npn4),
1967  & ti21(npn6,npn4,npn4),ti22(npn6,npn4,npn4)
1968  COMPLEX*16 CIM(NPN1)
1969 
1970  COMMON /TMAT/ TR11,TR12,TR21,TR22,TI11,TI12,TI21,TI22
1971  COMMON /CBESS/ B1R,B1I,B2R,B2I
1972  COMMON /SS/ SSIGN
1973  equivalence( plus1(1),tr11(1,1,1) )
1974  equivalence(d1(1,1,1),plus1(1)),
1975  & (d2(1,1,1),plus1(npl1*npn4*npn4+1)),
1976  & (d3(1,1,1),plus1(npl1*npn4*npn4*2+1)),
1977  & (d4(1,1,1),plus1(npl1*npn4*npn4*3+1)),
1978  & (d5r(1,1,1),plus1(npl1*npn4*npn4*4+1))
1979 
1980  CALL fact
1981  CALL signum
1982  lmax=2*nmax
1983  l1max=lmax+1
1984  ci=(0d0,1d0)
1985  cim(1)=ci
1986  DO 2 i=2,nmax
1987  cim(i)=cim(i-1)*ci
1988  2 CONTINUE
1989  ssi(1)=1d0
1990  DO 3 i=1,lmax
1991  i1=i+1
1992  si=dfloat(2*i+1)
1993  ssi(i1)=si
1994  IF(i.LE.nmax) ssj(i)=dsqrt(si)
1995  3 CONTINUE
1996  ci=-ci
1997  DO 5 i=1,nmax
1998  si=ssj(i)
1999  cci=cim(i)
2000  DO 4 j=1,nmax
2001  sj=1d0/ssj(j)
2002  ccj=cim(j)*sj/cci
2003  fr(j,i)=ccj
2004  fi(j,i)=ccj*ci
2005  ff(j,i)=si*sj
2006  4 CONTINUE
2007  5 CONTINUE
2008  nmax1=nmax+1
2009 
2010 C ***** CALCULATION OF THE ARRAYS B1 AND B2 *****
2011 
2012  k1=1
2013  k2=0
2014  k3=0
2015  k4=1
2016  k5=1
2017  k6=2
2018 
2019 C PRINT 3300, B1,B2
2020  3300 FORMAT (' B1 AND B2')
2021  DO 100 n=1,nmax
2022 
2023 C ***** CALCULATION OF THE ARRAYS T1 AND T2 *****
2024 
2025 
2026  DO 10 nn=1,nmax
2027  m1max=min0(n,nn)+1
2028  DO 6 m1=1,m1max
2029  m=m1-1
2030  l1=npn6+m
2031  tt1=tr11(m1,n,nn)
2032  tt2=tr12(m1,n,nn)
2033  tt3=tr21(m1,n,nn)
2034  tt4=tr22(m1,n,nn)
2035  tt5=ti11(m1,n,nn)
2036  tt6=ti12(m1,n,nn)
2037  tt7=ti21(m1,n,nn)
2038  tt8=ti22(m1,n,nn)
2039  t1=tt1+tt2
2040  t2=tt3+tt4
2041  t3=tt5+tt6
2042  t4=tt7+tt8
2043  tr1(l1,nn)=t1+t2
2044  tr2(l1,nn)=t1-t2
2045  ti1(l1,nn)=t3+t4
2046  ti2(l1,nn)=t3-t4
2047  IF(m.EQ.0) GO TO 6
2048  l1=npn6-m
2049  t1=tt1-tt2
2050  t2=tt3-tt4
2051  t3=tt5-tt6
2052  t4=tt7-tt8
2053  tr1(l1,nn)=t1-t2
2054  tr2(l1,nn)=t1+t2
2055  ti1(l1,nn)=t3-t4
2056  ti2(l1,nn)=t3+t4
2057  6 CONTINUE
2058  10 CONTINUE
2059 
2060 C ***** END OF THE CALCULATION OF THE ARRAYS T1 AND T2 *****
2061 
2062  nn1max=nmax1+n
2063  DO 40 nn1=1,nn1max
2064  n1=nn1-1
2065 
2066 C ***** CALCULATION OF THE ARRAYS A1 AND A2 *****
2067 
2068  CALL ccg(n,n1,nmax,k1,k2,g1)
2069  nnmax=min0(nmax,n1+n)
2070  nnmin=max0(1,iabs(n-n1))
2071  kn=n+nn1
2072  DO 15 nn=nnmin,nnmax
2073  nnn=nn+1
2074  sig=ssign(kn+nn)
2075  m1max=min0(n,nn)+npn6
2076  aar1=0d0
2077  aar2=0d0
2078  aai1=0d0
2079  aai2=0d0
2080  DO 13 m1=npn6,m1max
2081  m=m1-npn6
2082  sss=g1(m1,nnn)
2083  rr1=tr1(m1,nn)
2084  ri1=ti1(m1,nn)
2085  rr2=tr2(m1,nn)
2086  ri2=ti2(m1,nn)
2087  IF(m.EQ.0) GO TO 12
2088  m2=npn6-m
2089  rr1=rr1+tr1(m2,nn)*sig
2090  ri1=ri1+ti1(m2,nn)*sig
2091  rr2=rr2+tr2(m2,nn)*sig
2092  ri2=ri2+ti2(m2,nn)*sig
2093  12 aar1=aar1+sss*rr1
2094  aai1=aai1+sss*ri1
2095  aar2=aar2+sss*rr2
2096  aai2=aai2+sss*ri2
2097  13 CONTINUE
2098  xr=fr(nn,n)
2099  xi=fi(nn,n)
2100  ar1(nn)=aar1*xr-aai1*xi
2101  ai1(nn)=aar1*xi+aai1*xr
2102  ar2(nn)=aar2*xr-aai2*xi
2103  ai2(nn)=aar2*xi+aai2*xr
2104  15 CONTINUE
2105 
2106 C ***** END OF THE CALCULATION OF THE ARRAYS A1 AND A2 ****
2107 
2108  CALL ccg(n,n1,nmax,k3,k4,g2)
2109  m1=max0(-n1+1,-n)
2110  m2=min0(n1+1,n)
2111  m1max=m2+npn6
2112  m1min=m1+npn6
2113  DO 30 m1=m1min,m1max
2114  bbr1=0d0
2115  bbi1=0d0
2116  bbr2=0d0
2117  bbi2=0d0
2118  DO 25 nn=nnmin,nnmax
2119  nnn=nn+1
2120  sss=g2(m1,nnn)
2121  bbr1=bbr1+sss*ar1(nn)
2122  bbi1=bbi1+sss*ai1(nn)
2123  bbr2=bbr2+sss*ar2(nn)
2124  bbi2=bbi2+sss*ai2(nn)
2125  25 CONTINUE
2126  b1r(nn1,m1,n)=bbr1
2127  b1i(nn1,m1,n)=bbi1
2128  b2r(nn1,m1,n)=bbr2
2129  b2i(nn1,m1,n)=bbi2
2130  30 CONTINUE
2131  40 CONTINUE
2132  100 CONTINUE
2133 
2134 C ***** END OF THE CALCULATION OF THE ARRAYS B1 AND B2 ****
2135 
2136 C ***** CALCULATION OF THE ARRAYS D1,D2,D3,D4, AND D5 *****
2137 
2138 C PRINT 3301
2139  3301 FORMAT(' D1, D2, ...')
2140  DO 200 n=1,nmax
2141  DO 190 nn=1,nmax
2142  m1=min0(n,nn)
2143  m1max=npn6+m1
2144  m1min=npn6-m1
2145  nn1max=nmax1+min0(n,nn)
2146  DO 180 m1=m1min,m1max
2147  m=m1-npn6
2148  nn1min=iabs(m-1)+1
2149  dd1=0d0
2150  dd2=0d0
2151  DO 150 nn1=nn1min,nn1max
2152  xx=ssi(nn1)
2153  x1=b1r(nn1,m1,n)
2154  x2=b1i(nn1,m1,n)
2155  x3=b1r(nn1,m1,nn)
2156  x4=b1i(nn1,m1,nn)
2157  x5=b2r(nn1,m1,n)
2158  x6=b2i(nn1,m1,n)
2159  x7=b2r(nn1,m1,nn)
2160  x8=b2i(nn1,m1,nn)
2161  dd1=dd1+xx*(x1*x3+x2*x4)
2162  dd2=dd2+xx*(x5*x7+x6*x8)
2163  150 CONTINUE
2164  d1(m1,nn,n)=dd1
2165  d2(m1,nn,n)=dd2
2166  180 CONTINUE
2167  mmax=min0(n,nn+2)
2168  mmin=max0(-n,-nn+2)
2169  m1max=npn6+mmax
2170  m1min=npn6+mmin
2171  DO 186 m1=m1min,m1max
2172  m=m1-npn6
2173  nn1min=iabs(m-1)+1
2174  dd3=0d0
2175  dd4=0d0
2176  dd5r=0d0
2177  dd5i=0d0
2178  m2=-m+2+npn6
2179  DO 183 nn1=nn1min,nn1max
2180  xx=ssi(nn1)
2181  x1=b1r(nn1,m1,n)
2182  x2=b1i(nn1,m1,n)
2183  x3=b2r(nn1,m1,n)
2184  x4=b2i(nn1,m1,n)
2185  x5=b1r(nn1,m2,nn)
2186  x6=b1i(nn1,m2,nn)
2187  x7=b2r(nn1,m2,nn)
2188  x8=b2i(nn1,m2,nn)
2189  dd3=dd3+xx*(x1*x5+x2*x6)
2190  dd4=dd4+xx*(x3*x7+x4*x8)
2191  dd5r=dd5r+xx*(x3*x5+x4*x6)
2192  dd5i=dd5i+xx*(x4*x5-x3*x6)
2193  183 CONTINUE
2194  d3(m1,nn,n)=dd3
2195  d4(m1,nn,n)=dd4
2196  d5r(m1,nn,n)=dd5r
2197  d5i(m1,nn,n)=dd5i
2198  186 CONTINUE
2199  190 CONTINUE
2200  200 CONTINUE
2201 
2202 C ***** END OF THE CALCULATION OF THE D-ARRAYS *****
2203 
2204 C ***** CALCULATION OF THE EXPANSION COEFFICIENTS *****
2205 
2206 C PRINT 3303
2207  3303 FORMAT (' G1, G2, ...')
2208 
2209  dk=lam*lam/(4d0*csca*dacos(-1d0))
2210  DO 300 l1=1,l1max
2211  g1l=0d0
2212  g2l=0d0
2213  g3l=0d0
2214  g4l=0d0
2215  g5lr=0d0
2216  g5li=0d0
2217  l=l1-1
2218  sl=ssi(l1)*dk
2219  DO 290 n=1,nmax
2220  nnmin=max0(1,iabs(n-l))
2221  nnmax=min0(nmax,n+l)
2222  IF(nnmax.LT.nnmin) GO TO 290
2223  CALL ccg(n,l,nmax,k1,k2,g1)
2224  IF(l.GE.2) CALL ccg(n,l,nmax,k5,k6,g2)
2225  nl=n+l
2226  DO 280 nn=nnmin,nnmax
2227  nnn=nn+1
2228  mmax=min0(n,nn)
2229  m1min=npn6-mmax
2230  m1max=npn6+mmax
2231  si=ssign(nl+nnn)
2232  dm1=0d0
2233  dm2=0d0
2234  DO 270 m1=m1min,m1max
2235  m=m1-npn6
2236  IF(m.GE.0) sss1=g1(m1,nnn)
2237  IF(m.LT.0) sss1=g1(npn6-m,nnn)*si
2238  dm1=dm1+sss1*d1(m1,nn,n)
2239  dm2=dm2+sss1*d2(m1,nn,n)
2240  270 CONTINUE
2241  ffn=ff(nn,n)
2242  sss=g1(npn6+1,nnn)*ffn
2243  g1l=g1l+sss*dm1
2244  g2l=g2l+sss*dm2*si
2245  IF(l.LT.2) GO TO 280
2246  dm3=0d0
2247  dm4=0d0
2248  dm5r=0d0
2249  dm5i=0d0
2250  mmax=min0(n,nn+2)
2251  mmin=max0(-n,-nn+2)
2252  m1max=npn6+mmax
2253  m1min=npn6+mmin
2254  DO 275 m1=m1min,m1max
2255  m=m1-npn6
2256  sss1=g2(npn6-m,nnn)
2257  dm3=dm3+sss1*d3(m1,nn,n)
2258  dm4=dm4+sss1*d4(m1,nn,n)
2259  dm5r=dm5r+sss1*d5r(m1,nn,n)
2260  dm5i=dm5i+sss1*d5i(m1,nn,n)
2261  275 CONTINUE
2262  g5lr=g5lr-sss*dm5r
2263  g5li=g5li-sss*dm5i
2264  sss=g2(npn4,nnn)*ffn
2265  g3l=g3l+sss*dm3
2266  g4l=g4l+sss*dm4*si
2267  280 CONTINUE
2268  290 CONTINUE
2269  g1l=g1l*sl
2270  g2l=g2l*sl
2271  g3l=g3l*sl
2272  g4l=g4l*sl
2273  g5lr=g5lr*sl
2274  g5li=g5li*sl
2275  alf1(l1)=g1l+g2l
2276  alf2(l1)=g3l+g4l
2277  alf3(l1)=g3l-g4l
2278  alf4(l1)=g1l-g2l
2279  bet1(l1)=g5lr*2d0
2280  bet2(l1)=g5li*2d0
2281  lmax=l
2282  IF(abs(g1l).LT.1d-6) GO TO 500
2283  300 CONTINUE
2284  500 CONTINUE
2285  RETURN
2286  END
2287 
2288 C****************************************************************
2289 
2290 C CALCULATION OF THE QUANTITIES F(N+1)=0.5*LN(N!)
2291 C 0.LE.N.LE.899
2292 
2293  SUBROUTINE fact
2294  real*8 f(900)
2295  COMMON /fac/ f
2296  f(1)=0d0
2297  f(2)=0d0
2298  DO 2 i=3,900
2299  i1=i-1
2300  f(i)=f(i1)+0.5d0*dlog(dfloat(i1))
2301  2 CONTINUE
2302  RETURN
2303  END
2304 
2305 C************************************************************
2306 
2307 C CALCULATION OF THE ARRAY SSIGN(N+1)=SIGN(N)
2308 C 0.LE.N.LE.899
2309 
2310  SUBROUTINE signum
2311  real*8 ssign(900)
2312  COMMON /ss/ ssign
2313  ssign(1)=1d0
2314  DO 2 n=2,899
2315  ssign(n)=-ssign(n-1)
2316  2 CONTINUE
2317  RETURN
2318  END
2319 
2320 C******************************************************************
2321 C
2322 C Calculation of Clebsch-Gordan coefficients
2323 C (n,m:n1,m1/nn,mm)
2324 C for given n and n1. m1=mm-m, index mm is found from m as
2325 C mm=m*k1+k2
2326 C
2327 C Input parameters : N,N1,NMAX,K1,K2
2328 C N.LE.NMAX
2329 C N.GE.1
2330 C N1.GE.0
2331 C N1.LE.N+NMAX
2332 C Output parameters : GG(M+NPN6,NN+1) - array of the corresponding
2333 C coefficients
2334 C /M/.LE.N
2335 C /M1/=/M*(K1-1)+K2/.LE.N1
2336 C NN.LE.MIN(N+N1,NMAX)
2337 C NN.GE.MAX(/MM/,/N-N1/)
2338 C If K1=1 and K2=0, then 0.LE.M.LE.N
2339 
2340  SUBROUTINE ccg(N,N1,NMAX,K1,K2,GG)
2341  include 'tm_tmd.par.f'
2342  IMPLICIT REAL*8 (a-h,o-z)
2343  real*8 gg(npl1,npn6),cd(0:npn5),cu(0:npn5)
2344  IF(nmax.LE.npn4.
2345  & and.0.LE.n1.
2346  & and.n1.LE.nmax+n.
2347  & and.n.GE.1.
2348  & and.n.LE.nmax) GO TO 1
2349  print 5001
2350  stop
2351  5001 FORMAT(' ERROR IN SUBROUTINE CCG')
2352  1 nnf=min0(n+n1,nmax)
2353  min=npn6-n
2354  mf=npn6+n
2355  IF(k1.EQ.1.AND.k2.EQ.0) min=npn6
2356  DO 100 mind=min,mf
2357  m=mind-npn6
2358  mm=m*k1+k2
2359  m1=mm-m
2360  IF(iabs(m1).GT.n1) GO TO 90
2361  nnl=max0(iabs(mm),iabs(n-n1))
2362  IF(nnl.GT.nnf) GO TO 90
2363  nnu=n+n1
2364  nnm=(nnu+nnl)*0.5d0
2365  IF (nnu.EQ.nnl) nnm=nnl
2366  CALL ccgin(n,n1,m,mm,c)
2367  cu(nnl)=c
2368  IF (nnl.EQ.nnf) GO TO 50
2369  c2=0d0
2370  c1=c
2371  DO 7 nn=nnl+1,min0(nnm,nnf)
2372  a=dfloat((nn+mm)*(nn-mm)*(n1-n+nn))
2373  a=a*dfloat((n-n1+nn)*(n+n1-nn+1)*(n+n1+nn+1))
2374  a=dfloat(4*nn*nn)/a
2375  a=a*dfloat((2*nn+1)*(2*nn-1))
2376  a=dsqrt(a)
2377  b=0.5d0*dfloat(m-m1)
2378  d=0d0
2379  IF(nn.EQ.1) GO TO 5
2380  b=dfloat(2*nn*(nn-1))
2381  b=dfloat((2*m-mm)*nn*(nn-1)-mm*n*(n+1)+
2382  & mm*n1*(n1+1))/b
2383  d=dfloat(4*(nn-1)*(nn-1))
2384  d=d*dfloat((2*nn-3)*(2*nn-1))
2385  d=dfloat((nn-mm-1)*(nn+mm-1)*(n1-n+nn-1))/d
2386  d=d*dfloat((n-n1+nn-1)*(n+n1-nn+2)*(n+n1+nn))
2387  d=dsqrt(d)
2388  5 c=a*(b*c1-d*c2)
2389  c2=c1
2390  c1=c
2391  cu(nn)=c
2392  7 CONTINUE
2393  IF (nnf.LE.nnm) GO TO 50
2394  CALL direct(n,m,n1,m1,nnu,mm,c)
2395  cd(nnu)=c
2396  IF (nnu.EQ.nnm+1) GO TO 50
2397  c2=0d0
2398  c1=c
2399  DO 12 nn=nnu-1,nnm+1,-1
2400  a=dfloat((nn-mm+1)*(nn+mm+1)*(n1-n+nn+1))
2401  a=a*dfloat((n-n1+nn+1)*(n+n1-nn)*(n+n1+nn+2))
2402  a=dfloat(4*(nn+1)*(nn+1))/a
2403  a=a*dfloat((2*nn+1)*(2*nn+3))
2404  a=dsqrt(a)
2405  b=dfloat(2*(nn+2)*(nn+1))
2406  b=dfloat((2*m-mm)*(nn+2)*(nn+1)-mm*n*(n+1)
2407  & +mm*n1*(n1+1))/b
2408  d=dfloat(4*(nn+2)*(nn+2))
2409  d=d*dfloat((2*nn+5)*(2*nn+3))
2410  d=dfloat((nn+mm+2)*(nn-mm+2)*(n1-n+nn+2))/d
2411  d=d*dfloat((n-n1+nn+2)*(n+n1-nn-1)*(n+n1+nn+3))
2412  d=dsqrt(d)
2413  c=a*(b*c1-d*c2)
2414  c2=c1
2415  c1=c
2416  cd(nn)=c
2417  12 CONTINUE
2418  50 DO 9 nn=nnl,nnf
2419  IF (nn.LE.nnm) gg(mind,nn+1)=cu(nn)
2420  IF (nn.GT.nnm) gg(mind,nn+1)=cd(nn)
2421 c WRITE (6,*) N,M,N1,M1,NN,MM,GG(MIND,NN+1)
2422  9 CONTINUE
2423  90 CONTINUE
2424  100 CONTINUE
2425  RETURN
2426  END
2427 
2428 C*********************************************************************
2429 
2430  SUBROUTINE direct (N,M,N1,M1,NN,MM,C)
2431  IMPLICIT REAL*8 (a-h,o-z)
2432  real*8 f(900)
2433  COMMON /fac/ f
2434  c=f(2*n+1)+f(2*n1+1)+f(n+n1+m+m1+1)+f(n+n1-m-m1+1)
2435  c=c-f(2*(n+n1)+1)-f(n+m+1)-f(n-m+1)-f(n1+m1+1)-f(n1-m1+1)
2436  c=dexp(c)
2437  RETURN
2438  END
2439 
2440 C*********************************************************************
2441 C
2442 C Calculation of the Clebcsh-Gordan coefficients
2443 C G=(n,m:n1,mm-m/nn,mm)
2444 C for given n,n1,m,mm, where NN=MAX(/MM/,/N-N1/)
2445 C /M/.LE.N
2446 C /MM-M/.LE.N1
2447 C /MM/.LE.N+N1
2448 
2449  SUBROUTINE ccgin(N,N1,M,MM,G)
2450  IMPLICIT REAL*8 (a-h,o-z)
2451  real*8 f(900),ssign(900)
2452  COMMON /ss/ ssign
2453  COMMON /fac/ f
2454  m1=mm-m
2455  IF(n.GE.iabs(m).
2456  & and.n1.GE.iabs(m1).
2457  & and.iabs(mm).LE.(n+n1)) GO TO 1
2458  print 5001
2459  stop
2460  5001 FORMAT(' ERROR IN SUBROUTINE CCGIN')
2461  1 IF (iabs(mm).GT.iabs(n-n1)) GO TO 100
2462  l1=n
2463  l2=n1
2464  l3=m
2465  IF(n1.LE.n) GO TO 50
2466  k=n
2467  n=n1
2468  n1=k
2469  k=m
2470  m=m1
2471  m1=k
2472  50 n2=n*2
2473  m2=m*2
2474  n12=n1*2
2475  m12=m1*2
2476  g=ssign(n1+m1+1)
2477  & *dexp(f(n+m+1)+f(n-m+1)+f(n12+1)+f(n2-n12+2)-f(n2+2)
2478  & -f(n1+m1+1)-f(n1-m1+1)-f(n-n1+mm+1)-f(n-n1-mm+1))
2479  n=l1
2480  n1=l2
2481  m=l3
2482  RETURN
2483  100 a=1d0
2484  l1=m
2485  l2=mm
2486  IF(mm.GE.0) GO TO 150
2487  mm=-mm
2488  m=-m
2489  m1=-m1
2490  a=ssign(mm+n+n1+1)
2491  150 g=a*ssign(n+m+1)
2492  & *dexp(f(2*mm+2)+f(n+n1-mm+1)+f(n+m+1)+f(n1+m1+1)
2493  & -f(n+n1+mm+2)-f(n-n1+mm+1)-f(-n+n1+mm+1)-f(n-m+1)
2494  & -f(n1-m1+1))
2495  m=l1
2496  mm=l2
2497  RETURN
2498  END
2499 
2500 C*****************************************************************
2501 
2502  SUBROUTINE sarea (D,RAT)
2503  IMPLICIT REAL*8 (a-h,o-z)
2504  IF (d.GE.1) GO TO 10
2505  e=dsqrt(1d0-d*d)
2506  r=0.5d0*(d**(2d0/3d0) + d**(-1d0/3d0)*dasin(e)/e)
2507  r=dsqrt(r)
2508  rat=1d0/r
2509  RETURN
2510  10 e=dsqrt(1d0-1d0/(d*d))
2511  r=0.25d0*(2d0*d**(2d0/3d0) + d**(-4d0/3d0)*dlog((1d0+e)/(1d0-e))
2512  & /e)
2513  r=dsqrt(r)
2514  rat=1d0/r
2515  RETURN
2516  END
2517 
2518 c****************************************************************
2519 
2520  SUBROUTINE surfch (N,E,RAT)
2521  IMPLICIT REAL*8 (a-h,o-z)
2522  real*8 x(60),w(60)
2523  dn=dfloat(n)
2524  e2=e*e
2525  en=e*dn
2526  ng=60
2527  CALL gauss (ng,0,0,x,w)
2528  s=0d0
2529  v=0d0
2530  DO 10 i=1,ng
2531  xi=x(i)
2532  dx=dacos(xi)
2533  dxn=dn*dx
2534  ds=dsin(dx)
2535  dsn=dsin(dxn)
2536  dcn=dcos(dxn)
2537  a=1d0+e*dcn
2538  a2=a*a
2539  ens=en*dsn
2540  s=s+w(i)*a*dsqrt(a2+ens*ens)
2541  v=v+w(i)*(ds*a+xi*ens)*ds*a2
2542  10 CONTINUE
2543  rs=dsqrt(s*0.5d0)
2544  rv=(v*3d0/4d0)**(1d0/3d0)
2545  rat=rv/rs
2546  RETURN
2547  END
2548 
2549 C********************************************************************
2550 
2551  SUBROUTINE sareac (EPS,RAT)
2552  IMPLICIT REAL*8 (a-h,o-z)
2553  rat=(1.5d0/eps)**(1d0/3d0)
2554  rat=rat/dsqrt( (eps+2d0)/(2d0*eps) )
2555  RETURN
2556  END
2557 
2558 C********************************************************************
2559 
2560 C Computation of R1 and R2 for a power law size distribution with
2561 C effective radius A and effective variance B
2562 
2563  SUBROUTINE power (A,B,R1,R2)
2564  IMPLICIT REAL*8 (a-h,o-z)
2565  EXTERNAL f
2566  COMMON aa,bb
2567  aa=a
2568  bb=b
2569  ax=1d-5
2570  bx=a-1d-5
2571  r1=zeroin(ax,bx,f,0d0)
2572  r2=(1d0+b)*2d0*a-r1
2573  RETURN
2574  END
2575 
2576 C***********************************************************************
2577 
2578  DOUBLE PRECISION FUNCTION zeroin (AX,BX,F,TOL)
2579  IMPLICIT REAL*8 (a-h,o-z)
2580  eps=1d0
2581  10 eps=0.5d0*eps
2582  tol1=1d0+eps
2583  IF (tol1.GT.1d0) GO TO 10
2584  15 a=ax
2585  b=bx
2586  fa=f(a)
2587  fb=f(b)
2588  20 c=a
2589  fc=fa
2590  d=b-a
2591  e=d
2592  30 IF (abs(fc).GE.abs(fb)) GO TO 40
2593  35 a=b
2594  b=c
2595  c=a
2596  fa=fb
2597  fb=fc
2598  fc=fa
2599  40 tol1=2d0*eps*abs(b)+0.5d0*tol
2600  xm=0.5d0*(c-b)
2601  IF (abs(xm).LE.tol1) GO TO 90
2602  44 IF (fb.EQ.0d0) GO TO 90
2603  45 IF (abs(e).LT.tol1) GO TO 70
2604  46 IF (abs(fa).LE.abs(fb)) GO TO 70
2605  47 IF (a.NE.c) GO TO 50
2606  48 s=fb/fa
2607  p=2d0*xm*s
2608  q=1d0-s
2609  GO TO 60
2610  50 q=fa/fc
2611  r=fb/fc
2612  s=fb/fa
2613  p=s*(2d0*xm*q*(q-r)-(b-a)*(r-1d0))
2614  q=(q-1d0)*(r-1d0)*(s-1d0)
2615  60 IF (p.GT.0d0) q=-q
2616  p=abs(p)
2617  IF ((2d0*p).GE.(3d0*xm*q-abs(tol1*q))) GO TO 70
2618  64 IF (p.GE.abs(0.5d0*e*q)) GO TO 70
2619  65 e=d
2620  d=p/q
2621  GO TO 80
2622  70 d=xm
2623  e=d
2624  80 a=b
2625  fa=fb
2626  IF (abs(d).GT.tol1) b=b+d
2627  IF (abs(d).LE.tol1) b=b+dsign(tol1,xm)
2628  fb=f(b)
2629  IF ((fb*(fc/abs(fc))).GT.0d0) GO TO 20
2630  85 GO TO 30
2631  90 zeroin=b
2632  RETURN
2633  END
2634 
2635 C***********************************************************************
2636 
2637  DOUBLE PRECISION FUNCTION f(R1)
2638  IMPLICIT REAL*8 (a-h,o-z)
2639  COMMON a,b
2640  r2=(1d0+b)*2d0*a-r1
2641  f=(r2-r1)/dlog(r2/r1)-a
2642  RETURN
2643  END
2644 
2645 C**********************************************************************
2646 C CALCULATION OF POINTS AND WEIGHTS OF GAUSSIAN QUADRATURE *
2647 C FORMULA. IF IND1 = 0 - ON INTERVAL (-1,1), IF IND1 = 1 - ON *
2648 C INTERVAL (0,1). IF IND2 = 1 RESULTS ARE PRINTED. *
2649 C N - NUMBER OF POINTS *
2650 C Z - DIVISION POINTS *
2651 C W - WEIGHTS *
2652 C**********************************************************************
2653 
2654  SUBROUTINE gauss ( N,IND1,IND2,Z,W )
2655  IMPLICIT REAL*8 (a-h,p-z)
2656  real*8 z(n),w(n)
2657  a=1d0
2658  b=2d0
2659  c=3d0
2660  ind=mod(n,2)
2661  k=n/2+ind
2662  f=dfloat(n)
2663  DO 100 i=1,k
2664  m=n+1-i
2665  IF(i.EQ.1) x=a-b/((f+a)*f)
2666  IF(i.EQ.2) x=(z(n)-a)*4d0+z(n)
2667  IF(i.EQ.3) x=(z(n-1)-z(n))*1.6d0+z(n-1)
2668  IF(i.GT.3) x=(z(m+1)-z(m+2))*c+z(m+3)
2669  IF(i.EQ.k.AND.ind.EQ.1) x=0d0
2670  niter=0
2671  check=1d-16
2672  10 pb=1d0
2673  niter=niter+1
2674  IF (niter.LE.100) GO TO 15
2675  check=check*10d0
2676  15 pc=x
2677  dj=a
2678  DO 20 j=2,n
2679  dj=dj+a
2680  pa=pb
2681  pb=pc
2682  20 pc=x*pb+(x*pb-pa)*(dj-a)/dj
2683  pa=a/((pb-x*pc)*f)
2684  pb=pa*pc*(a-x*x)
2685  x=x-pb
2686  IF(dabs(pb).GT.check*dabs(x)) GO TO 10
2687  z(m)=x
2688  w(m)=pa*pa*(a-x*x)
2689  IF(ind1.EQ.0) w(m)=b*w(m)
2690  IF(i.EQ.k.AND.ind.EQ.1) GO TO 100
2691  z(i)=-z(m)
2692  w(i)=w(m)
2693  100 CONTINUE
2694  IF(ind2.NE.1) GO TO 110
2695  print 1100,n
2696  1100 FORMAT(' *** POINTS AND WEIGHTS OF GAUSSIAN QUADRATURE FORMULA',
2697  * ' OF ',i4,'-TH ORDER')
2698  DO 105 i=1,k
2699  zz=-z(i)
2700  105 print 1200,i,zz,i,w(i)
2701  1200 FORMAT(' ',4x,'X(',i4,') = ',f17.14,5x,'W(',i4,') = ',f17.14)
2702  GO TO 115
2703  110 CONTINUE
2704 C PRINT 1300,N
2705  1300 FORMAT(' GAUSSIAN QUADRATURE FORMULA OF ',i4,'-TH ORDER IS USED')
2706  115 CONTINUE
2707  IF(ind1.EQ.0) GO TO 140
2708  DO 120 i=1,n
2709  120 z(i)=(a+z(i))/b
2710  140 CONTINUE
2711  RETURN
2712  END
2713 
2714 C****************************************************************
2715 
2716  SUBROUTINE distrb (NNK,YY,WY,NDISTR,AA,BB,GAM,R1,R2,REFF,
2717  & VEFF,PI)
2718  IMPLICIT REAL*8 (a-h,o-z)
2719  real*8 yy(nnk),wy(nnk)
2720  IF (ndistr.EQ.2) GO TO 100
2721  IF (ndistr.EQ.3) GO TO 200
2722  IF (ndistr.EQ.4) GO TO 300
2723  IF (ndistr.EQ.5) GO TO 360
2724  print 1001,aa,bb,gam
2725  1001 FORMAT('MODIFIED GAMMA DISTRIBUTION, alpha=',f6.4,' r_c=',
2726  & f6.4,' gamma=',f6.4)
2727  a2=aa/gam
2728  db=1d0/bb
2729  DO 50 i=1,nnk
2730  x=yy(i)
2731  y=x**aa
2732  x=x*db
2733  y=y*dexp(-a2*(x**gam))
2734  wy(i)=wy(i)*y
2735  50 CONTINUE
2736  GO TO 400
2737  100 print 1002,aa,bb
2738  1002 FORMAT('LOG-NORMAL DISTRIBUTION, r_g=',f8.4,
2739  & ' [ln(sigma_g)]**2=', f6.4)
2740  da=1d0/aa
2741  DO 150 i=1,nnk
2742  x=yy(i)
2743  y=dlog(x*da)
2744  y=dexp(-y*y*0.5d0/bb)/x
2745  wy(i)=wy(i)*y
2746  150 CONTINUE
2747  GO TO 400
2748  200 print 1003
2749  1003 FORMAT('POWER LAW DISTRIBUTION OF HANSEN & TRAVIS 1974')
2750  DO 250 i=1,nnk
2751  x=yy(i)
2752  wy(i)=wy(i)/(x*x*x)
2753  250 CONTINUE
2754  GO TO 400
2755  300 print 1004,aa,bb
2756  1004 FORMAT ('GAMMA DISTRIBUTION, a=',f6.3,' b=',f6.4)
2757  b2=(1d0-3d0*bb)/bb
2758  dab=1d0/(aa*bb)
2759  DO 350 i=1,nnk
2760  x=yy(i)
2761  x=(x**b2)*dexp(-x*dab)
2762  wy(i)=wy(i)*x
2763  350 CONTINUE
2764  GO TO 400
2765  360 print 1005,bb
2766  1005 FORMAT ('MODIFIED POWER LAW DISTRIBUTION, alpha=',d10.4)
2767  DO 370 i=1,nnk
2768  x=yy(i)
2769  IF (x.LE.r1) wy(i)=wy(i)
2770  IF (x.GT.r1) wy(i)=wy(i)*(x/r1)**bb
2771  370 CONTINUE
2772  400 CONTINUE
2773  sum=0d0
2774  DO 450 i=1,nnk
2775  sum=sum+wy(i)
2776  450 CONTINUE
2777  sum=1d0/sum
2778  DO 500 i=1,nnk
2779  wy(i)=wy(i)*sum
2780  500 CONTINUE
2781  g=0d0
2782  DO 550 i=1,nnk
2783  x=yy(i)
2784  g=g+x*x*wy(i)
2785  550 CONTINUE
2786  reff=0d0
2787  DO 600 i=1,nnk
2788  x=yy(i)
2789  reff=reff+x*x*x*wy(i)
2790  600 CONTINUE
2791  reff=reff/g
2792  veff=0d0
2793  DO 650 i=1,nnk
2794  x=yy(i)
2795  xi=x-reff
2796  veff=veff+xi*xi*x*x*wy(i)
2797  650 CONTINUE
2798  veff=veff/(g*reff*reff)
2799  RETURN
2800  END
2801 
2802 C*************************************************************
2803 
2804  SUBROUTINE hovenr(L1,A1,A2,A3,A4,B1,B2)
2805  IMPLICIT REAL*8 (a-h,o-z)
2806  real*8 a1(l1),a2(l1),a3(l1),a4(l1),b1(l1),b2(l1)
2807  DO 100 l=1,l1
2808  kontr=1
2809  ll=l-1
2810  dl=dfloat(ll)*2d0+1d0
2811  ddl=dl*0.48d0
2812  aa1=a1(l)
2813  aa2=a2(l)
2814  aa3=a3(l)
2815  aa4=a4(l)
2816  bb1=b1(l)
2817  bb2=b2(l)
2818  IF(ll.GE.1.AND.abs(aa1).GE.dl) kontr=2
2819  IF(abs(aa2).GE.dl) kontr=2
2820  IF(abs(aa3).GE.dl) kontr=2
2821  IF(abs(aa4).GE.dl) kontr=2
2822  IF(abs(bb1).GE.ddl) kontr=2
2823  IF(abs(bb2).GE.ddl) kontr=2
2824  IF(kontr.EQ.2) print 3000,ll
2825  c=-0.1d0
2826  DO 50 i=1,11
2827  c=c+0.1d0
2828  cc=c*c
2829  c1=cc*bb2*bb2
2830  c2=c*aa4
2831  c3=c*aa3
2832  IF((dl-c*aa1)*(dl-c*aa2)-cc*bb1*bb1.LE.-1d-4) kontr=2
2833  IF((dl-c2)*(dl-c3)+c1.LE.-1d-4) kontr=2
2834  IF((dl+c2)*(dl-c3)-c1.LE.-1d-4) kontr=2
2835  IF((dl-c2)*(dl+c3)-c1.LE.-1d-4) kontr=2
2836  IF(kontr.EQ.2) print 4000,ll,c
2837  50 CONTINUE
2838  100 CONTINUE
2839 C IF(KONTR.EQ.1) PRINT 2000
2840  2000 FORMAT('TEST OF VAN DER MEE & HOVENIER IS SATISFIED')
2841  3000 FORMAT('TEST OF VAN DER MEE & HOVENIER IS NOT SATISFIED, L=',i3)
2842  4000 FORMAT('TEST OF VAN DER MEE & HOVENIER IS NOT SATISFIED, L=',i3,
2843  & ' A=',d9.2)
2844  RETURN
2845  END
2846 
2847 C****************************************************************
2848 
2849 C CALCULATION OF THE SCATTERING MATRIX FOR GIVEN EXPANSION
2850 C COEFFICIENTS
2851 
2852 C A1,...,B2 - EXPANSION COEFFICIENTS
2853 C LMAX - NUMBER OF COEFFICIENTS MINUS 1
2854 C N - NUMBER OF SCATTERING ANGLES
2855 C THE CORRESPONDING SCATTERING ANGLES ARE GIVEN BY
2856 C 180*(I-1)/(N-1) (DEGREES), WHERE I NUMBERS THE ANGLES
2857 
2858  SUBROUTINE matr(A1,A2,A3,A4,B1,B2,LMAX,NPNA,NPNAX,IAX,FM)
2859  include 'tm_tmd.par.f'
2860  IMPLICIT REAL*8 (a-h,o-z)
2861  real*8 a1(npl),a2(npl),a3(npl),a4(npl),b1(npl),b2(npl),fm(6,npna,npnax)
2862  n=npna
2863  dn=1d0/dfloat(n-1)
2864  da=dacos(-1d0)*dn
2865  db=180d0*dn
2866  l1max=lmax+1
2867 C PRINT 1000
2868  1000 FORMAT(' ')
2869 C PRINT 1001
2870  1001 FORMAT(' ',2x,'S',6x,'ALPHA1',6x,'ALPHA2',6x,'ALPHA3',
2871  & 6x,'ALPHA4',7x,'BETA1',7x,'BETA2')
2872  DO 10 l1=1,l1max
2873  l=l1-1
2874 C PRINT 1002,L,A1(L1),A2(L1),A3(L1),A4(L1),B1(L1),B2(L1)
2875  10 CONTINUE
2876  1002 FORMAT(' ',i3,6f12.5)
2877  tb=-db
2878  taa=-da
2879  print 1000
2880  print 1003
2881  1003 FORMAT(' ',5x,'<',8x,'F11',8x,'F22',8x,'F33',
2882  & 8x,'F44',8x,'F12',8x,'F34')
2883  d6=dsqrt(6d0)*0.25d0
2884  DO 500 i1=1,n
2885  taa=taa+da
2886  tb=tb+db
2887  u=dcos(taa)
2888  f11=0d0
2889  f2=0d0
2890  f3=0d0
2891  f44=0d0
2892  f12=0d0
2893  f34=0d0
2894  p1=0d0
2895  p2=0d0
2896  p3=0d0
2897  p4=0d0
2898  pp1=1d0
2899  pp2=0.25d0*(1d0+u)*(1d0+u)
2900  pp3=0.25d0*(1d0-u)*(1d0-u)
2901  pp4=d6*(u*u-1d0)
2902  DO 400 l1=1,l1max
2903  l=l1-1
2904  dl=dfloat(l)
2905  dl1=dfloat(l1)
2906  f11=f11+a1(l1)*pp1
2907  f44=f44+a4(l1)*pp1
2908  IF(l.EQ.lmax) GO TO 350
2909  pl1=dfloat(2*l+1)
2910  p=(pl1*u*pp1-dl*p1)/dl1
2911  p1=pp1
2912  pp1=p
2913  350 IF(l.LT.2) GO TO 400
2914  f2=f2+(a2(l1)+a3(l1))*pp2
2915  f3=f3+(a2(l1)-a3(l1))*pp3
2916  f12=f12+b1(l1)*pp4
2917  f34=f34+b2(l1)*pp4
2918  IF(l.EQ.lmax) GO TO 400
2919  pl2=dfloat(l*l1)*u
2920  pl3=dfloat(l1*(l*l-4))
2921  pl4=1d0/dfloat(l*(l1*l1-4))
2922  p=(pl1*(pl2-4d0)*pp2-pl3*p2)*pl4
2923  p2=pp2
2924  pp2=p
2925  p=(pl1*(pl2+4d0)*pp3-pl3*p3)*pl4
2926  p3=pp3
2927  pp3=p
2928  p=(pl1*u*pp4-dsqrt(dfloat(l*l-4))*p4)/dsqrt(dfloat(l1*l1-4))
2929  p4=pp4
2930  pp4=p
2931  400 CONTINUE
2932  f22=(f2+f3)*0.5d0
2933  f33=(f2-f3)*0.5d0
2934 C F22=F22/F11
2935 C F33=F33/F11
2936 C F44=F44/F11
2937 C F12=-F12/F11
2938 C F34=F34/F11
2939  print 1004,tb,f11,f22,f33,f44,f12,f34
2940  fm(1,i1,iax)=f11
2941  fm(2,i1,iax)=f22
2942  fm(3,i1,iax)=f33
2943  fm(4,i1,iax)=f44
2944  fm(5,i1,iax)=f12
2945  fm(6,i1,iax)=f34
2946  500 CONTINUE
2947  print 1000
2948  1004 FORMAT(' ',f6.2,6f11.4)
2949  RETURN
2950  END
vary
subroutine vary(LAM, MRR, MRI, A, EPS, NP, NGAUSS, X, P, PPI, PIR, PII, R, DR, DDR, DRR, DRI, NMAX)
Definition: ampld.lp.f:978
f1
float f1(float x)
rsp1
subroutine rsp1(X, NG, NGAUSS, REV, EPS, NP, R, DR)
Definition: ampld.lp.f:1030
zeroin
double precision function zeroin(AX, BX, F, TOL)
Definition: tmd.lp.f:1395
hovenr
subroutine hovenr(L1, A1, A2, A3, A4, B1, B2)
Definition: tmd.lp.f:1561
ccg
subroutine ccg(N, N1, NMAX, K1, K2, GG)
Definition: tmd.lp.f:1209
tmatr0
subroutine tmatr0(NGAUSS, X, W, AN, ANN, S, SS, PPI, PIR, PII, R, DR, DDR, DRR, DRI, NMAX, NCHECK)
Definition: ampld.lp.f:1301
rjb
subroutine rjb(X, Y, U, NMAX, NNMAX)
Definition: ampld.lp.f:1172
gsp
subroutine gsp(NMAX, CSCA, LAM, ALF1, ALF2, ALF3, ALF4, BET1, BET2, LMAX)
Definition: tmd.lp.f:813
f3
float f3(float z)
tmatr
subroutine tmatr(M, NGAUSS, X, W, AN, ANN, S, SS, PPI, PIR, PII, R, DR, DDR, DRR, DRI, NMAX, NCHECK)
Definition: ampld.lp.f:1514
real
#define real
Definition: DbAlgOcean.cpp:26
const
subroutine const(NGAUSS, NMAX, MMAX, P, X, W, AN, ANN, S, SS, NP, EPS)
Definition: ampld.lp.f:924
rsp3
subroutine rsp3(X, NG, NGAUSS, REV, EPS, R, DR)
Definition: ampld.lp.f:1078
signum
subroutine signum
Definition: tmd.lp.f:1178
vig
subroutine vig(X, NMAX, M, DV1, DV2)
Definition: ampld.lp.f:1793
fac
#define fac
Definition: atmocor1_land.c:111
f2
float f2(float y)
matr
subroutine matr(A1, A2, A3, A4, B1, B2, LMAX, NPNA, F)
Definition: tmd.lp.f:1615
sareac
subroutine sareac(EPS, RAT)
Definition: ampld.lp.f:1951
sarea
subroutine sarea(D, RAT)
Definition: ampld.lp.f:1902
surfch
subroutine surfch(N, E, RAT)
Definition: ampld.lp.f:1920
zgetri
subroutine zgetri(N, A, LDA, IPIV, WORK, LWORK, INFO)
Definition: lpd.f:1160
zgetrf
subroutine zgetrf(M, N, A, LDA, IPIV, INFO)
Definition: lpd.f:138
max
#define max(A, B)
Definition: main_biosmap.c:61
tt
subroutine tt(NMAX, NCHECK)
Definition: ampld.lp.f:1852
and
===========================================================================V5.0.48(Terra) 03/20/2015 Changes shown below are differences from MOD_PR02 V5.0.46(Terra)============================================================================Changes noted for V6.1.20(Terra) below were also instituted for this version.============================================================================V6.1.20(Terra) 03/12/2015 Changes shown below are differences from MOD_PR02 V6.1.18(Terra)============================================================================Changes from v6.1.18 which may affect scientific output:A situation can occur in which a scan which contains sector rotated data has a telemetry value indicating the completeness of the sector rotation. This issue is caused by the timing of the instrument command to perform the sector rotation and the recording of the telemetry point that reports the status of sector rotation. In this case a scan is considered valid by L1B and pass through the calibration - reporting extremely high radiances. Operationally the TEB calibration uses a 40 scan average coefficient, so the 20 scans(one mirror side) after the sector rotation are contaminated with anomalously high radiance values. A similar timing issue appeared before the sector rotation was fixed in V6.1.2. Our analysis indicates the ‘SET_FR_ENC_DELTA’ telemetry correlates well with the sector rotation encoder position. The use of this telemetry point to determine scans that are sector rotated should fix the anomaly occured before and after the sector rotation(usually due to the lunar roll maneuver). The fix related to the sector rotation in V6.1.2 is removed in this version.============================================================================V6.1.18(Terra) 10/01/2014 Changes shown below are differences from MOD_PR02 V6.1.16(Terra)============================================================================Added doi attributes to NRT(Near-Real-Time) product.============================================================================V6.1.16(Terra) 01/27/2014 Changes shown below are differences from MOD_PR02 V6.1.14(Terra)============================================================================Migrate to SDP Toolkit 5.2.17============================================================================V6.1.14(Terra) 06/26/2012 Changes shown below are differences from MOD_PR02 V6.1.12(Terra)============================================================================Added the doi metadata to L1B product============================================================================V6.1.12(Terra) 04/25/2011 Changes shown below are differences from MOD_PR02 V6.1.8(Terra)============================================================================1. The algorithm to calculate uncertainties for reflective solar bands(RSB) is updated. The current uncertainty in L1B code includes 9 terms from prelaunch analysis. The new algorithm regroups them with the new added contributions into 5 terms:u1:the common term(AOI and time independent) and
Definition: HISTORY.txt:126
fact
subroutine fact
Definition: tmd.lp.f:1161
min
#define min(A, B)
Definition: main_biosmap.c:62
rsp2
subroutine rsp2(X, NG, REV, EPS, N, R, DR)
Definition: ampld.lp.f:1053
pi
#define pi
Definition: vincenty.c:23
bess
subroutine bess(X, XR, XI, NG, NMAX, NNMAX1, NNMAX2)
Definition: ampld.lp.f:1138
ryb
subroutine ryb(X, Y, V, NMAX)
Definition: ampld.lp.f:1199
power
subroutine power(A, B, R1, R2)
Definition: tmd.lp.f:1380
ccgin
subroutine ccgin(N, N1, M, MM, G)
Definition: tmd.lp.f:1318
calcrand
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)
Definition: tmd.lp.f:362
distrb
subroutine distrb(NNK, YY, WY, NDISTR, AA, BB, GAM, R1, R2, REFF, VEFF, PI)
Definition: tmd.lp.f:1474
f
#define f
Definition: l1_czcs_hdf.c:702
l3
Definition: L3File.cpp:10
gauss
subroutine gauss(x1, x2, x, w, n)
Definition: 6sm1.f:4055
abs
#define abs(a)
Definition: misc.h:90
direct
subroutine direct(N, M, N1, M1, NN, MM, C)
Definition: tmd.lp.f:1299
cjb
subroutine cjb(XR, XI, YR, YI, UR, UI, NMAX, NNMAX)
Definition: ampld.lp.f:1228
der
subroutine der(T, X, DX)
Definition: der.f:2


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Mission Quality Monitor

HELP

Level 1 & 2 Browser Help

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility Policy

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster