2 this result was generated by the
code in its current setting.
5 randomly oriented oblate spheroids, a/b= 2.0000000
6 lam= .500000 mrr= .1330d+01 mri= .1000d-02
7 accuracy of computations ddelt = .10d-02
8 power law distribution of hansen & travis 1974
9 r1= .546765 r2= 1.653235
10 equal-surface-area-sphere reff= 1.0000 veff= .1000
11 number of gaussian quadrature points in
SIZE averaging = 7
12 test of van
der mee & hovenier is satisfied
13 cext= .620059d+01 csca= .608811d+01 w= .981861d+00 <cos>= .795853d+00
15 s alpha1 alpha2 alpha3 alpha4 beta1 beta2
16 0 1.00000 .00000 .00000 .94701 .00000 .00000
17 1 2.38756 .00000 .00000 2.40317 .00000 .00000
18 2 3.44047 4.22550 4.14121 3.44978 .00518 .14199
19 3 4.30784 4.52938 4.47654 4.30075 -.06700 .00061
20 4 4.80175 5.22540 5.18611 4.77897 -.07461 -.04857
21 5 5.01955 5.30004 5.29356 5.01550 -.03166 -.00095
22 6 5.20044 5.40802 5.40301 5.21102 -.04963 .03442
23 7 5.38064 5.54781 5.52966 5.38320 -.07723 -.03423
24 8 5.44501 5.61122 5.58415 5.43366 -.06418 -.07369
25 9 5.38994 5.53497 5.52494 5.39340 -.04911 -.06812
26 10 5.30885 5.43390 5.42481 5.31525 -.06113 -.07484
27 11 5.17481 5.29188 5.27444 5.17434 -.06655 -.11320
28 12 4.97178 5.08303 5.06276 4.96487 -.04977 -.12667
29 13 4.72347 4.81759 4.81172 4.72934 -.03833 -.12035
30 14 4.46612 4.55484 4.54481 4.46935 -.04693 -.12703
31 15 4.18891 4.26669 4.25410 4.18946 -.04134 -.14121
32 16 3.89903 3.97221 3.95599 3.89420 -.03223 -.14110
33 17 3.60608 3.66490 3.65715 3.60738 -.01908 -.12943
34 18 3.33839 3.39259 3.38341 3.33836 -.01783 -.12134
35 19 3.07906 3.12521 3.11950 3.08222 -.01682 -.11384
36 20 2.83593 2.87894 2.86939 2.83473 -.02121 -.11290
37 21 2.59089 2.62667 2.61851 2.59005 -.01798 -.10778
38 22 2.35901 2.39386 2.38536 2.35686 -.01420 -.10238
39 23 2.13152 2.15936 2.15669 2.13463 -.01221 -.09111
40 24 1.92719 1.95493 1.94712 1.92542 -.01460 -.09125
41 25 1.72730 1.74913 1.74417 1.72716 -.00846 -.08019
42 26 1.55255 1.57218 1.56664 1.55139 -.01261 -.07069
43 27 1.38960 1.40551 1.40243 1.39042 -.01655 -.06411
44 28 1.24209 1.25881 1.25242 1.23950 -.01708 -.06463
45 29 1.09356 1.10609 1.10493 1.09581 -.01618 -.05606
46 30 .96469 .97814 .97202 .96161 -.01519 -.05859
47 31 .83777 .84762 .84634 .83921 -.00963 -.04780
48 32 .73508 .74432 .74017 .73340 -.01009 -.04482
49 33 .63980 .64648 .64486 .64002 -.00898 -.03741
50 34 .55901 .56592 .56303 .55828 -.01224 -.03509
51 35 .48236 .48737 .48754 .48451 -.01467 -.03293
52 36 .41721 .42356 .41848 .41401 -.01304 -.03843
53 37 .34674 .35079 .35055 .34808 -.00841 -.02956
54 38 .29376 .29824 .29524 .29220 -.00766 -.03002
55 39 .24175 .24475 .24412 .24220 -.00422 -.02486
56 40 .19831 .20150 .19936 .19725 -.00225 -.02311
57 41 .15699 .15882 .15929 .15815 -.00105 -.01718
58 42 .12641 .12892 .12638 .12466 -.00144 -.01831
59 43 .09392 .09522 .09581 .09518 -.00006 -.01285
60 44 .07194 .07377 .07139 .07019 -.00016 -.01455
61 45 .04664 .04752 .04759 .04714 .00131 -.00887
62 46 .03170 .03279 .03082 .03018 .00028 -.00796
63 47 .01815 .01858 .01833 .01809 .00085 -.00418
64 48 .01168 .01213 .01091 .01067 .00087 -.00397
65 49 .00612 .00629 .00548 .00539 .00137 -.00204
66 50 .00362 .00375 .00254 .00246 .00117 -.00132
67 51 .00132 .00137 .00135 .00131 .00042 -.00029
68 52 .00080 .00084 .00072 .00070 .00035 -.00025
69 53 .00041 .00042 .00029 .00028 .00028 -.00012
70 54 .00023 .00023 .00010 .00009 .00020 -.00005
71 55 -.00001 -.00001 -.00002 -.00002 .00000 .00003
72 56 .00000 .00000 .00000 .00000 .00000 .00000
73 57 .00000 .00000 .00000 .00000 .00000 .00000
74 58 .00000 .00000 .00000 .00000 .00000 .00000
76 < f11 f22 f33 f44 f12 f34
77 .00 110.2658 110.2341 110.2341 110.2025 .0000 .0000
78 10.00 23.6825 23.6703 23.6349 23.6304 .2337 .7009
79 20.00 4.8125 4.8052 4.7649 4.7664 .1925 .0423
80 30.00 2.0213 2.0137 1.9813 1.9836 .0881 -.0595
81 40.00 .9749 .9679 .9405 .9429 .0558 -.0385
82 50.00 .4846 .4782 .4519 .4540 .0340 -.0372
83 60.00 .2596 .2534 .2239 .2257 .0327 -.0370
84 70.00 .2006 .1931 .1535 .1550 .0132 -.0644
85 80.00 .2069 .1978 .1490 .1496 .0072 -.0952
86 90.00 .2490 .2375 .1704 .1714 -.0266 -.1384
87 100.00 .2368 .2217 .1662 .1692 -.0474 -.1096
88 110.00 .1793 .1595 .1168 .1251 -.0316 -.0589
89 120.00 .1412 .1126 .0653 .0832 -.0163 -.0415
90 130.00 .1201 .0788 .0298 .0600 -.0059 -.0309
91 140.00 .1045 .0532 .0041 .0447 .0008 -.0229
92 150.00 .0882 .0383 -.0142 .0281 .0008 -.0175
93 160.00 .0791 .0348 -.0208 .0189 .0002 -.0128
94 170.00 .0754 .0468 -.0171 .0080 .0056 -.0061
95 180.00 .1404 .0653 -.0653 .0097 .0000 .0000
97 power law distribution of hansen & travis 1974
98 r1= .273383 r2= .826617
99 equal-surface-area-sphere reff= .5000 veff= .1000
100 number of gaussian quadrature points in
SIZE averaging = 4
101 test of van
der mee & hovenier is satisfied
102 cext= .173748d+01 csca= .172372d+01 w= .992080d+00 <cos>= .837254d+00
104 s alpha1 alpha2 alpha3 alpha4 beta1 beta2
105 0 1.00000 .00000 .00000 .96133 .00000 .00000
106 1 2.51176 .00000 .00000 2.51910 .00000 .00000
107 2 3.49636 4.33696 4.27436 3.48836 .02198 .06483
108 3 4.06342 4.52880 4.50112 4.05338 -.06414 -.01285
109 4 4.23806 4.71668 4.69771 4.23276 -.05219 -.07709
110 5 4.13236 4.50467 4.49117 4.13400 -.03610 -.08912
111 6 3.86520 4.17435 4.15638 3.86421 -.05188 -.10175
112 7 3.50844 3.75244 3.73735 3.50750 -.04031 -.12583
113 8 3.11901 3.31220 3.30076 3.12128 -.03132 -.11254
114 9 2.72759 2.87978 2.87048 2.73212 -.04365 -.11080
115 10 2.34961 2.47532 2.46112 2.34827 -.04691 -.12001
116 11 1.97824 2.07695 2.06801 1.98063 -.02986 -.11621
117 12 1.64863 1.72767 1.71815 1.64894 -.02985 -.09926
118 13 1.35245 1.41252 1.40499 1.35315 -.02988 -.09096
119 14 1.08940 1.13961 1.13177 1.08836 -.02730 -.08448
120 15 .84868 .88749 .88297 .84915 -.02075 -.07558
121 16 .63643 .66904 .66340 .63513 -.01352 -.06865
122 17 .45239 .47584 .47284 .45262 -.00289 -.05384
123 18 .31306 .33021 .32625 .31133 -.00111 -.04100
124 19 .20413 .21430 .21251 .20383 .00108 -.02688
125 20 .13219 .14030 .13779 .13084 -.00050 -.02051
126 21 .07694 .08129 .08041 .07647 .00141 -.01300
127 22 .04280 .04582 .04471 .04212 .00096 -.00843
128 23 .02029 .02194 .02162 .02022 .00148 -.00402
129 24 .01036 .01130 .01084 .01002 .00049 -.00274
130 25 .00285 .00311 .00299 .00276 .00072 -.00085
131 26 .00086 .00095 .00090 .00081 .00026 -.00023
132 27 .00020 .00023 .00021 .00019 .00008 -.00004
133 28 .00003 .00004 .00003 .00003 .00002 .00000
134 29 .00000 .00000 .00000 .00000 .00000 .00000
135 30 .00000 .00000 .00000 .00000 .00000 .00000
137 < f11 f22 f33 f44 f12 f34
138 .00 43.8217 43.7953 43.7953 43.7689 .0000 .0000
139 10.00 25.7096 25.6858 25.6793 25.6622 .1209 .4479
140 20.00 7.6985 7.6819 7.6544 7.6474 .1490 .3276
141 30.00 2.7576 2.7466 2.7192 2.7169 .0704 .1134
142 40.00 1.0872 1.0790 1.0549 1.0541 .0577 .0262
143 50.00 .5419 .5354 .5122 .5122 .0322 -.0192
144 60.00 .3062 .3003 .2766 .2770 .0200 -.0436
145 70.00 .2236 .2175 .1937 .1938 -.0045 -.0575
146 80.00 .1727 .1657 .1409 .1414 -.0147 -.0499
147 90.00 .1361 .1285 .0999 .1009 -.0182 -.0523
148 100.00 .1147 .1048 .0702 .0730 -.0315 -.0525
149 110.00 .0982 .0855 .0481 .0536 -.0446 -.0360
150 120.00 .0799 .0634 .0270 .0366 -.0368 -.0143
151 130.00 .0643 .0438 .0058 .0206 -.0206 -.0065
152 140.00 .0574 .0345 -.0067 .0109 -.0128 -.0055
153 150.00 .0480 .0287 -.0138 .0013 -.0057 -.0051
154 160.00 .0432 .0289 -.0152 -.0036 -.0009 -.0073
155 170.00 .0404 .0277 -.0211 -.0100 .0049 .0009
156 180.00 .0623 .0318 -.0318 -.0013 .0000 .0000
160 c new release including the lapack
matrix inversion procedure.
161 c we thank cory davis(university of edinburgh) for pointing
162 c out the possibility of replacing the proprietary nag
matrix
163 c inversion routine by the
public-domain lapack equivalent.
165 c calculation of light scattering by polydisperse, randomly
166 c oriented particles of identical axially symmetric shape
168 c this version of the
code uses extended precision variables
169 c
and must be used along with the accompanying files tmq.par.
f
172 c last update 08/06/2005
174 c the
code has been developed by michael mishchenko at the nasa
175 c goddard institute for space studies, new york. this research
176 c was funded by the nasa radiation sciences program.
178 c the
code can be used without limitations in any not-for-
179 c profit scientific research. we only request that in any
180 c publication using the
code the source of the
code be acknowledged
181 c
and relevant references(see below) be made.
183 c this version of the
code is applicable to spheroids,
184 c chebyshev particles,
and finite circular cylinders.
186 c the computational method is based
on the watermsn
's T-matrix
187 C approach and is described in detail in the following papers:
189 C 1. M. I. Mishchenko, Light scattering by randomly oriented
190 C axially symmetric particles, J. Opt. Soc. Am. A,
191 C vol. 8, 871-882 (1991).
193 C 2. M. I. Mishchenko, Light scattering by size-shape
194 C distributions of randomly oriented axially symmetric
195 C particles of a size comparable to a wavelength,
196 C Appl. Opt., vol. 32, 4652-4666 (1993).
198 C 3. M. I. Mishchenko and L. D. Travis, T-matrix computations
199 C of light scattering by large spheroidal particles,
200 C Opt. Commun., vol. 109, 16-21 (1994).
202 C 4. M. I. Mishchenko, L. D. Travis, and A. Macke, Scattering
203 C of light by polydisperse, randomly oriented, finite
204 C circular cylinders, Appl. Opt., vol. 35, 4927-4940 (1996).
206 C 5. D. J. Wielaard, M. I. Mishchenko, A. Macke, and B. E. Carlson,
207 C Improved T-matrix computations for large, nonabsorbing and
208 C weakly absorbing nonspherical particles and comparison
209 C with geometrical optics approximation, Appl. Opt., vol. 36,
212 C A general review of the T-matrix approach can be found in
214 C 6. M. I. Mishchenko, L. D. Travis, and D. W. Mackowski,
215 C T-matrix computations of light scattering by nonspherical
216 C particles: a review, J. Quant. Spectrosc. Radiat.
217 C Transfer, vol. 55, 535-575 (1996).
219 C The following paper provides a detailed user guide to the
222 C 7. M. I. Mishchenko and L. D. Travis, Capabilities and
223 C limitations of a current FORTRAN implementation of the
224 C T-matrix method for randomly oriented, rotationally
225 C symmetric scatterers, J. Quant. Spectrosc. Radiat. Transfer,
226 C vol. 60, 309-324 (1998).
228 C These papers are available in the .pdf format at the web site
230 C http://www.giss.nasa.gov/~crmim/publications/
232 C or in hardcopy upon request from Michael Mishchenko
233 C Please e-mail your request to crmim@giss.nasa.gov.
235 C A comprehensive book "Scattering, Absorption, and Emission of
236 C Light by Small Particles" (Cambridge University Press, Cambridge,
237 C 2002) is also available in the .pdf format at the web site
239 C http://www.giss.nasa.gov/~crmim/books.html
241 C Analytical averaging over particle orientations (Ref. 1) makes
242 C this method the fastest exact technique currently available.
243 C The use of an automatic convergence procedure
244 C (Ref. 2) makes the code convenient in massive computations.
245 C Ref. 4 describes features specific for finite cylinders as
246 C particles with sharp rectangular edges. Ref. 5 describes further
247 C numerical improvements.
249 C Since this code uses extended precision variables, it is slower
250 C than the equivalent double-precision code. The double-precision
251 C code is also available, but it cannot compute as large particles
252 C for the same shape, refractive index, and wavelength as this code.
254 C This is the first part of the full T-matrix code. The second part,
255 C lpq.f, is completely independent of the first part. It contains no
256 C T-matrix-specific subroutines and can be compiled separately.
257 C The second part of the code replaces the previously implemented
258 C standard matrix inversion scheme based on Gaussian elimination
259 C by a scheme based on the LU factorization technique.
260 C As described in Ref. 5 above, the use of the LU factorization is
261 C especially beneficial for nonabsorbing or weakly absorbing particles.
262 C In this code we use the LAPACK implementation of the LU factorization
263 C scheme. LAPACK stands for Linear Algebra PACKage. The latter is
264 C publicly available at the following internet site:
266 C http://www.netlib.org/lapack/
270 C RAT = 1 - particle size is specified in terms of the
271 C equal-volume-sphere radius
272 .NE.
C RAT1 - particle size is specified in terms of the
273 C equal-surface-area-sphere radius
274 C NDISTR specifies the distribution of equivalent-sphere radii
275 C NDISTR = 1 - modified gamma distribution
276 C [Eq. (40) of Ref. 7]
280 C NDISTR = 2 - log-normal distribution
281 C [Eq. 41) of Ref. 7]
284 C NDISTR = 3 - power law distribution
285 C [Eq. (42) of Ref. 7]
286 C AXI=r_eff (effective radius)
287 C B=v_eff (effective variance)
288 C Parameters R1 and R2 (see below) are calculated
289 C automatically for given AXI and B
290 C NDISTR = 4 - gamma distribution
291 C [Eq. (39) of Ref. 7]
294 C NDISTR = 5 - modified power law distribution
295 C [Eq. (24) in M. I. Mishchenko et al.,
296 C Bidirectional reflectance of flat,
297 C optically thick particulate laters: an efficient radiative
298 C transfer solution and applications to snow and soil surfaces,
299 C J. Quant. Spectrosc. Radiat. Transfer, Vol. 63, 409-432 (1999)].
302 C The code computes NPNAX size distributions of the same type
303 C and with the same values of B and GAM in one run.
304 C The parameter AXI varies from AXMAX to AXMAX/NPNAX in steps of
305 C AXMAX/NPNAX. To compute a single size distribution, use
306 C NPNAX=1 and AXMAX equal to AXI of this size distribution.
308 C R1 and R2 - minimum and maximum equivalent-sphere
309 C radii in the size distribution.
310 C They are calculated automatically
311 C for the power law distribution with given AXI and B
312 C but must be specified for other distributions
316 C AXI=AXMAX-DAX*DFLOAT(IAX-1)
318 C in the main program.
319 C For the modified power law distribution (NDISTR=5), the
320 C minimum radius is 0, R2 is the maximum radius,
321 C and R1 is the intermediate radius at which the
322 C n(r)=const dependence is replaced by the power law
325 .LE.
C NKMAX988 is such that NKMAX+2 is the
326 C number of Gaussian quadrature points used in
327 C integrating over the size distribution for particles with
328 C AXI=AXMAX. For particles with AXI=AXMAX-AXMAX/NPNAX,
329 C AXMAX-2*AXMAX/NPNAX, etc. the number of Gaussian points
330 C linearly decreases.
331 C For the modified power law distribution, the number
332 C of integration points on the interval [0,R1] is also
335 C LAM - wavelength of light
336 C MRR and MRI - real and imaginary parts of the refractive
338 C EPS and NP - specify the shape of the particles.
339 C For spheroids NP=-1 and EPS is the ratio of the
340 C horizontal to rotational axes. EPS is larger than
341 C 1 for oblate spheroids and smaller than 1 for
343 C For cylinders NP=-2 and EPS is the ratio of the
344 C diameter to the length.
345 C For Chebyshev particles NP must be positive and
346 C is the degree of the Chebyshev polynomial, while
347 C EPS is the deformation parameter
348 C [Eq. (33) of Ref. 7].
349 C DDELT - accuracy of the computations
350 C NPNA - number of equidistant scattering angles (from 0
351 C to 180 deg) for which the scattering matrix is
353 C NDGS - parameter controlling the number of division points
354 C in computing integrals over the particle surface.
355 C For compact particles, the recommended value is 2.
356 C For highly aspherical particles larger values (3, 4,...)
357 C may be necessary to obtain convergence.
358 C The code does not check convergence over this parameter.
359 C Therefore, control comparisons of results obtained with
360 C different NDGS-values are recommended.
365 C REFF and VEFF - effective radius and effective variance of
366 C the size distribution as defined by Eqs. (43)-(45) of
368 C CEXT - extinction cross section per particle
369 C CSCA - scattering cross section per particle
370 C W - single scattering albedo
371 C <cos> - asymmetry parameter of the phase function
372 C ALPHA1,...,BETA2 - coefficients appearing in the expansions
373 C of the elements of the scattering matrix in
374 C generalized spherical functions
375 C [Eqs. (11)-(16) of Ref. 7].
376 C F11,...,F44 - elements of the normalized scattering matrix [as
377 C defined by Eqs. (1)-(3) of Ref. 7] versus scattering angle
379 C Note that LAM, r_c, r_g, r_eff, a, R1, and R2 must
380 C be given in the same units of length, and that
381 C the dimension of CEXT and CSCA is that of LAM squared (e.g., if
382 C LAM and AXI are given in microns, then CEXT and CSCA are
383 C calculated in square microns).
385 C The physical correctness of the computed results is tested using
386 C the general inequalities derived by van der Mee and Hovenier,
387 C Astron. Astrophys., vol. 228, 559-568 (1990). Although
388 C the message that the test of van der Mee and Hovenier is satisfied
389 C does not guarantee that the results are absolutely correct,
390 C the message that the test is not satisfied can mean that something
393 C The convergence of the T-matrix method for particles with
394 C different sizes, refractive indices, and aspect ratios can be
395 C dramatically different. Usually, large sizes and large aspect
396 C ratios cause problems. The user of this code
397 C should first experiment with different input parameters in
398 C order to get an idea of the range of applicability of this
399 C technique. Sometimes decreasing the aspect ratio
400 C from 3 to 2 can increase the maximum convergent equivalent-
401 C sphere size parameter by a factor of several (Ref. 7).
402 C The CPU time required rapidly increases with increasing ratio
403 C radius/wavelength and/or with increasing particle asphericity.
404 C This should be taken into account in planning massive computations.
405 C Using an optimizing compiler on IBM RISC workstations saves
406 C about 70% of CPU time.
408 C Execution can be automatically terminated if dimensions of certain
409 C arrays are not big enough or if the convergence procedure decides
410 C that the accuracy of extended precision variables is insufficient
411 C to obtain a converged T-matrix solution for given particles.
412 C In all cases, a message appears explaining the cause of termination.
415 C "WARNING: W IS GREATER THAN 1"
416 C means that the single-scattering albedo exceeds the maximum
417 C possible value 1. If W is greater than 1 by more than
418 C DDELT, this message can be an indication of numerical
419 C instability caused by extreme values of particle parameters.
421 C The message "WARNING: NGAUSS=NPNG1" means that convergence over
422 C the parameter NG (see Ref. 2) cannot be obtained for the NPNG1
423 C value specified in the PARAMETER statement in the file tmq.par.f.
424 C Often this is not a serious problem, especially for compact
427 C Larger and/or more aspherical particles may require larger
428 C values of the parameters NPN1, NPN4, and NPNG1 in the file
429 C tmq.par.f. It is recommended to keep NPN1=NPN4+25 and
430 C NPNG1=3*NPN1. Note that the memory requirement increases
431 C as the third power of NPN4. If the memory of
432 C a computer is too small to accomodate the code in its current
433 C setting, the parameters NPN1, NPN4, and NPNG1 should be
434 C decreased. However, this will decrease the maximum size parameter
435 C that can be handled by the code.
437 C In some cases any increases of NPN1 will not make the T-matrix
438 C computations convergent. This means that the particle is just
439 C too "bad" (extreme size parameter and/or extreme aspect ratio
440 C and/or extreme refractive index; see Ref. 7).
441 C The main program contains several PRINT statements which are
442 C currently commentd out. If uncommented, these statements will
443 C produce numbers which show the convergence rate and can be
444 C used to determine whether T-matrix computations for given particle
445 C parameters will converge at all.
447 C Some of the common blocks are used to save memory rather than
448 C to transfer data. Therefore, if a compiler produces a warning
449 C message that the lengths of a common block are different in
450 C different subroutines, this is not a real problem.
452 C The recommended value of DDELT is 0.001. For bigger values,
453 C false convergence can be obtained.
455 C In computations for spheres use EPS=1.000001 instead of EPS=1.
456 C The use of EPS=1 can cause overflows in some rare cases.
458 C To calculate a monodisperse particle, use the options
466 C AXI=AXMAX-DAX*DFLOAT(IAX-1)
470 C where R is the equivalent-sphere radius.
472 C It is recommended to use the power law rather than the
473 C gamma size distribution, because in this case convergent solution
474 C can be obtained for larger REFF and VEFF assuming the same
475 C maximal R2 (Mishchenko and Travis, Appl. Opt., vol. 33, 7206-7225,
478 C For some compilers, DACOS must be raplaced by DARCOS and DASIN
481 C If many different size distributions are computed and the
482 C refractive index is fixed, then another approach can be more
483 C efficient than running this code many times. Specifically,
484 C scattering results should be computed for monodisperse particles
485 C with sizes ranging from essentially zero to some maximum value
486 C with a small step size (say, 0.02 microns). These results
487 C should be stored on disk and can be used along with spline
488 C interpolation to compute scattering by particles with intermediate
489 C sizes. Scattering patterns for monodisperse nonspherical
490 C particles in random orientation are (much) smoother than for
491 C monodisperse spheres, and spline interpolation usually gives good
492 C results. In this way, averaging over any size distribution is a
493 C matter of seconds. For more on size averaging, see Refs. 2 and 4.
495 C We would highly appreciate informing me of any problems encountered
496 C with this code. Please send your message to the following
497 C e-mail address: CRMIM@GISS.NASA.GOV.
499 C WHILE THE COMPUTER PROGRAM HAS BEEN TESTED FOR A VARIETY OF CASES,
500 C IT IS NOT INCONCEIVABLE THAT IT CONTAINS UNDETECTED ERRORS. ALSO,
501 C INPUT PARAMETERS CAN BE USED WHICH ARE OUTSIDE THE ENVELOPE OF
502 C VALUES FOR WHICH RESULTS ARE COMPUTED ACCURATELY. FOR THIS REASON,
503 C THE AUTHORS AND THEIR ORGANIZATION DISCLAIM ALL LIABILITY FOR
504 C ANY DAMAGES THAT MAY RESULT FROM THE USE OF THE PROGRAM.
507 IMPLICIT REAL*8 (A-H,O-Z)
509 REAL*16 LAM,MRR,MRI,X(NPNG2),W(NPNG2),S(NPNG2),SS(NPNG2),
510 * AN(NPN1),R(NPNG2),DR(NPNG2),PPI,PIR,PII,P,EPS,A,
511 * DDR(NPNG2),DRR(NPNG2),DRI(NPNG2),ANN(NPN1,NPN1)
512 REAL*8 XG(1000),WG(1000),TR1(NPN2,NPN2),TI1(NPN2,NPN2),
513 & ALPH1(NPL),ALPH2(NPL),ALPH3(NPL),ALPH4(NPL),BET1(NPL),
514 & BET2(NPL),XG1(2000),WG1(2000),
515 & AL1(NPL),AL2(NPL),AL3(NPL),AL4(NPL),BE1(NPL),BE2(NPL)
517 & RT11(NPN6,NPN4,NPN4),RT12(NPN6,NPN4,NPN4),
518 & RT21(NPN6,NPN4,NPN4),RT22(NPN6,NPN4,NPN4),
519 & IT11(NPN6,NPN4,NPN4),IT12(NPN6,NPN4,NPN4),
520 & IT21(NPN6,NPN4,NPN4),IT22(NPN6,NPN4,NPN4)
523 COMMON /TMAT/ RT11,RT12,RT21,RT22,IT11,IT12,IT21,IT22
527 C OPEN FILES *******************************************************
530 OPEN (10,FILE='tmatr.write
')
532 C INPUT DATA ********************************************************
551 .EQ..OR..EQ.
IF (NP-1NP-2) NCHECK=1
552 .GT..AND..EQ.
IF (NP0(-1)**NP1) NCHECK=1
553 WRITE (6,5454) NCHECK
554 5454 FORMAT ('ncheck=
',I1)
560 .GT..AND..EQ.
IF (DABS(RAT-1D0)1D-8NP-1) CALL SAREA (DEPS,RAT)
561 .GT..AND..GE.
IF (DABS(RAT-1D0)1D-8NP0) CALL SURFCH(NP,DEPS,RAT)
562 .GT..AND..EQ.
IF (DABS(RAT-1D0)1D-8NP-2) CALL SAREAC (DEPS,RAT)
564 8000 FORMAT ('rat=
',F8.6)
565 .EQ..AND..GE.
IF(NP-1EPS1D0) PRINT 7000,EPS
566 .EQ..AND..LT.
IF(NP-1EPS1D0) PRINT 7001,EPS
567 .GE.
IF(NP0) PRINT 7100,NP,EPS
568 .EQ..AND..GE.
IF(NP-2EPS1D0) PRINT 7150,EPS
569 .EQ..AND..LT.
IF(NP-2EPS1D0) PRINT 7151,EPS
570 PRINT 7400, LAM,MRR,MRI
572 7000 FORMAT('randomly oriented oblate spheroids, a/b=
',F11.7)
573 7001 FORMAT('randomly oriented prolate spheroids, a/b=
',F11.7)
574 7100 FORMAT('randomly oriented chebyshev particles, t
',
576 7150 FORMAT('randomly oriented oblate cylinders, d/l=
',F11.7)
577 7151 FORMAT('randomly oriented prolate cylinders, d/l=
',F11.7)
578 7200 FORMAT ('accuracy of computations ddelt =
',d8.2)
579 7400 FORMAT('lam=
',F10.6,3X,'mrr=
',D10.4,3X,'mri=
',D10.4)
582 AXI=AXMAX-DAX*DFLOAT(IAX-1)
585 NK=IDINT(AXI*NKMAX/AXMAX+2)
586 .GT.
IF (NK1000) PRINT 8001,NK
587 .GT.
IF (NK1000) STOP
588 .EQ.
IF (NDISTR3) CALL POWER (AXI,B,R1,R2)
589 8001 FORMAT ('nk=
',I4,' i.e., is greater than 1000.
',
590 & 'execution terminated.
')
591 CALL GAUSS (NK,0,0,XG,WG)
595 .EQ.
IF (NDISTR5) GO TO 3
611 4 CALL DISTRB (NK,XG1,WG1,NDISTR,AXI,B,GAM,R1,R2,
614 8002 FORMAT('r1=
',F10.6,' r2=
',F10.6)
615 .LE.
IF (DABS(RAT-1D0)1D-6) PRINT 8003, REFF,VEFF
616 .GT.
IF (DABS(RAT-1D0)1D-6) PRINT 8004, REFF,VEFF
617 8003 FORMAT('equal-volume-sphere reff=
',F8.4,' veff=
',F7.4)
618 8004 FORMAT('equal-surface-area-sphere reff=
',F8.4,
621 7250 FORMAT('number of gaussian quadrature points
',
622 & 'in
SIZE averaging =
',I4)
638 IXXX=XEV+4.05D0*XEV**0.333333D0
640 .GE.
IF (INM1NPN1) PRINT 7333, NPN1
641 .GE.
IF (INM1NPN1) STOP
642 7333 FORMAT('convergence is not obtained for npn1=
',I3,
643 & '. execution terminated
')
650 .GT.
IF (NGAUSSNPNG1) PRINT 7340, NGAUSS
651 .GT.
IF (NGAUSSNPNG1) STOP
652 7340 FORMAT('ngauss =
',I3,' i.e. is greater than npng1.
',
653 & ' execution terminated
')
654 7334 FORMAT(' nmax =
', I3,' dc2=
',D8.2,' dc1=
',D8.2)
655 7335 FORMAT(' nmax1 =
', I3,' dc2=
',D8.2,
657 CALL CONST(NGAUSS,NMAX,MMAX,P,X,W,AN,ANN,S,SS,NP,EPS)
658 CALL VARY(LAM,MRR,MRI,A,EPS,NP,NGAUSS,X,P,PPI,PIR,PII,R,
659 & DR,DDR,DRR,DRI,NMAX)
660 CALL TMATR0 (NGAUSS,X,W,AN,ANN,S,SS,PPI,PIR,PII,R,DR,
661 & DDR,DRR,DRI,NMAX,NCHECK)
671 QSCA=QSCA+DN1*(TR1NN*TR1NN+TI1NN*TI1NN
672 & +TR1NN1*TR1NN1+TI1NN1*TI1NN1)
673 QEXT=QEXT+(TR1NN+TR1NN1)*DN1
675 DSCA=DABS((QSCA1-QSCA)/QSCA)
676 DEXT=DABS((QEXT1-QEXT)/QEXT)
677 C PRINT 7334, NMAX,DSCA,DEXT
680 NMIN=DFLOAT(NMAX)/2D0+1D0
688 DQSCA=DN1*(TR1NN*TR1NN+TI1NN*TI1NN
689 & +TR1NN1*TR1NN1+TI1NN1*TI1NN1)
690 DQEXT=(TR1NN+TR1NN1)*DN1
691 DQSCA=DABS(DQSCA/QSCA)
692 DQEXT=DABS(DQEXT/QEXT)
694 .LE..AND..LE.
IF (DQSCADDELTDQEXTDDELT) GO TO 12
697 c PRINT 7335, NMAX1,DQSCA,DQEXT
698 .LE..AND..LE.
IF(DSCADDELTDEXTDDELT) GO TO 55
699 .EQ.
IF (NMANPN1) PRINT 7333, NPN1
700 .EQ.
IF (NMANPN1) STOP
703 .EQ.
IF (NGAUSSNPNG1) PRINT 7336
705 DO 150 NGAUS=NNNGGG,NPNG1
708 7336 FORMAT('warning: ngauss=npng1
')
709 7337 FORMAT(' ng=
',I3,' dc2=
',D8.2,' dc1=
',D8.2)
710 CALL CONST(NGAUSS,NMAX,MMAX,P,X,W,AN,ANN,S,SS,NP,EPS)
711 CALL VARY(LAM,MRR,MRI,A,EPS,NP,NGAUSS,X,P,PPI,PIR,PII,R,
712 & DR,DDR,DRR,DRI,NMAX)
713 CALL TMATR0 (NGAUSS,X,W,AN,ANN,S,SS,PPI,PIR,PII,R,DR,
714 & DDR,DRR,DRI,NMAX,NCHECK)
724 QSCA=QSCA+DN1*(TR1NN*TR1NN+TI1NN*TI1NN
725 & +TR1NN1*TR1NN1+TI1NN1*TI1NN1)
726 QEXT=QEXT+(TR1NN+TR1NN1)*DN1
728 DSCA=DABS((QSCA1-QSCA)/QSCA)
729 DEXT=DABS((QEXT1-QEXT)/QEXT)
730 c PRINT 7337, NGGG,DSCA,DEXT
733 .LE..AND..LE.
IF(DSCADDELTDEXTDDELT) GO TO 155
734 .EQ.
IF (NGAUSNPNG1) PRINT 7336
743 .GT.
IF (NMAX1NPN4) PRINT 7550, NMAX1
744 7550 FORMAT ('nmax1 =
',I3, ', i.e. greater than npn4.
',
745 & ' execution terminated
')
746 .GT.
IF (NMAX1NPN4) STOP
767 QSCA=QSCA+ZZ1*ZZ1+ZZ2*ZZ2+ZZ3*ZZ3+ZZ4*ZZ4
768 & +ZZ5*ZZ5+ZZ6*ZZ6+ZZ7*ZZ7+ZZ8*ZZ8
770 C PRINT 7800,0,DABS(QEXT),QSCA,NMAX
772 CALL TMATR(M,NGAUSS,X,W,AN,ANN,S,SS,PPI,PIR,PII,R,DR,
773 & DDR,DRR,DRI,NMAX,NCHECK)
800 QSC=QSC+(ZZ1*ZZ1+ZZ2*ZZ2+ZZ3*ZZ3+ZZ4*ZZ4
801 & +ZZ5*ZZ5+ZZ6*ZZ6+ZZ7*ZZ7+ZZ8*ZZ8)*2D0
810 C PRINT 7800,M,DABS(QXT),QSC,NMAX
811 7800 FORMAT(' m=
',I3,' qxt=
',D12.6,' qsc=
',D12.6,
814 COEFF1=LAM*LAM*0.5D0/P
817 c PRINT 7880, NMAX,NMAX1
818 7880 FORMAT ('nmax=
',I3,' nmax1=
',I3)
819 CALL GSP (NMAX1,CSCA,LAM,AL1,AL2,AL3,AL4,BE1,BE2,LMAX)
825 ALPH1(L1)=ALPH1(L1)+AL1(L1)*WGI
826 ALPH2(L1)=ALPH2(L1)+AL2(L1)*WGI
827 ALPH3(L1)=ALPH3(L1)+AL3(L1)*WGI
828 ALPH4(L1)=ALPH4(L1)+AL4(L1)*WGI
829 BET1(L1)=BET1(L1)+BE1(L1)*WGI
830 BET2(L1)=BET2(L1)+BE2(L1)*WGI
833 CEXTIN=CEXTIN+CEXT*WGII
834 C PRINT 6070, I,NMAX,NMAX1,NGAUSS
838 ALPH1(L1)=ALPH1(L1)/CSCAT
839 ALPH2(L1)=ALPH2(L1)/CSCAT
840 ALPH3(L1)=ALPH3(L1)/CSCAT
841 ALPH4(L1)=ALPH4(L1)/CSCAT
842 BET1(L1)=BET1(L1)/CSCAT
843 BET2(L1)=BET2(L1)/CSCAT
846 CALL HOVENR(L1MAX,ALPH1,ALPH2,ALPH3,ALPH4,BET1,BET2)
848 PRINT 9100,CEXTIN,CSCAT,WALB,ASYMM
849 9100 FORMAT('cext=
',D12.6,2X,'csca=
',D12.6,2X,
850 & 2X,'w=
',D12.6,2X,'<cos>=
',D12.6)
851 .GT.
IF (WALB1D0) PRINT 9111
852 9111 FORMAT ('warning: w is greater than 1
')
853 WRITE (10,580) WALB,L1MAX
855 WRITE (10,575) ALPH1(L),ALPH2(L),ALPH3(L),ALPH4(L),
861 CALL MATR (ALPH1,ALPH2,ALPH3,ALPH4,BET1,BET2,LMAX,NPNA)
864 TIME=DFLOAT(ITIME)/6000D0
866 1001 FORMAT (' time =
',F8.2,' min')
870 C**********************************************************************
872 C INPUT PARAMETERS: *
874 C NG = 2*NGAUSS - number of quadrature points on the *
875 .LE.
C interval (-1,1). NGAUSSNPNG1 *
876 .LE.
C NMAX,MMAX - maximum dimensions of the arrays. NMAXNPN1 *
880 C OUTPUT PARAMETERS: *
882 C X,W - points and weights of the quadrature formula *
884 C ANN(N1,N2) = (1/2)*sqrt((2*n1+1)*(2*n2+1)/(n1*(n1+1)*n2*(n2+1))) *
885 C S(I)=1/sin(arccos(x(i))) *
888 C**********************************************************************
890 SUBROUTINE CONST (NGAUSS,NMAX,MMAX,P,X,W,AN,ANN,S,SS,NP,EPS)
891 IMPLICIT REAL*16 (A-H,O-Z)
893 REAL*16 X(NPNG2),W(NPNG2),X1(NPNG1),W1(NPNG1),
894 * X2(NPNG1),W2(NPNG1),
895 * S(NPNG2),SS(NPNG2),
896 * AN(NPN1),ANN(NPN1,NPN1),DD(NPN1)
901 D=QSQRT(QFLOAT(2*N+1)/QFLOAT(NN))
909 .EQ.
IF (NP-2) GO TO 11
910 CALL QGAUSS(NG,0,0,X,W)
912 11 NG1=DFLOAT(NGAUSS)/2D0
915 CALL QGAUSS(NG1,0,0,X1,W1)
916 CALL QGAUSS(NG2,0,0,X2,W2)
918 W(I)=0.5Q0*(XX+1Q0)*W1(I)
919 X(I)=0.5Q0*(XX+1Q0)*X1(I)+0.5Q0*(XX-1Q0)
922 W(I+NG1)=-0.5Q0*XX*W2(I)
923 X(I+NG1)=-0.5Q0*XX*X2(I)+0.5Q0*XX
941 C***************************************************************
943 SUBROUTINE QGAUSS ( N,IND1,IND2,Z,W )
944 IMPLICIT REAL*16 (A-H,P-Z)
954 .EQ.
IF(I1) X=A-B/((F+A)*F)
955 .EQ.
IF(I2) X=(Z(N)-A)*4Q0+Z(N)
956 .EQ.
IF(I3) X=(Z(N-1)-Z(N))*1.6Q0+Z(N-1)
957 .GT.
IF(I3) X=(Z(M+1)-Z(M+2))*C+Z(M+3)
958 .EQ..AND..EQ.
IF(IKIND1) X=0Q0
963 .LE.
IF (NITER100) GO TO 15
972 20 PC=X*PB+(X*PB-PA)*(DJ-A)/DJ
976 .GT.
IF(QABS(PB)CHECK*QABS(X)) GO TO 10
979 .EQ.
IF(IND10) W(M)=B*W(M)
980 .EQ..AND..EQ.
IF(IKIND1) GO TO 100
984 5000 format ('qgauss does not converge, check=
',Q10.3)
985 .NE.
IF(IND21) GO TO 110
987 1100 FORMAT(' *** points
and weights of gaussian quadrature formula
',
988 * ' of
',I4,'-th order
')
991 105 PRINT 1200,I,ZZ,I,W(I)
992 1200 FORMAT(' ',4X,'x(
',I4,') =
',F17.14,5X,'w(
',I4,') =
',F17.14)
996 1300 FORMAT(' gaussian quadrature formula of
',I4,'-th order is used
')
998 .EQ.
IF(IND10) GO TO 140
1005 C**********************************************************************
1007 C INPUT PARAMETERS: *
1009 C LAM - wavelength of light *
1010 C MRR and MRI - real and imaginary parts of the refractive index *
1011 C A,EPS,NP - specify shape of the particle *
1012 C (see subroutines RSP1, RSP2, and RSP3) *
1013 C NG = NGAUSS*2 - number of gaussian quadrature points on the *
1015 C X - gaussian division points *
1018 C OUTPUT INFORMATION: *
1020 C PPI = PI**2 , where PI = (2*P)/LAM (wavenumber) *
1023 C R and DR - see subroutines RSP1, RSP2, and RSP3 *
1024 C DDR=1/(PI*SQRT(R)) *
1025 C DRR+I*DRI=DDR/(MRR+I*MRI) *
1026 C NMAX - dimension of T(m)-matrices *
1027 C arrays J,Y,JR,JI,DJ,DY,DJR,DJI are transferred through *
1028 C COMMON /CBESS/ - see subroutine BESS *
1030 C**********************************************************************
1032 SUBROUTINE VARY (LAM,MRR,MRI,A,EPS,NP,NGAUSS,X,P,PPI,PIR,PII,
1033 * R,DR,DDR,DRR,DRI,NMAX)
1035 IMPLICIT REAL*16 (A-H,O-Z)
1036 REAL*16 X(NPNG2),R(NPNG2),DR(NPNG2),MRR,MRI,LAM,
1037 * Z(NPNG2),ZR(NPNG2),ZI(NPNG2),
1038 * J(NPNG2,NPN1),Y(NPNG2,NPN1),JR(NPNG2,NPN1),
1039 * JI(NPNG2,NPN1),DJ(NPNG2,NPN1),
1040 * DJR(NPNG2,NPN1),DJI(NPNG2,NPN1),DDR(NPNG2),
1041 * DRR(NPNG2),DRI(NPNG2),
1043 COMMON /CBESS/ J,Y,JR,JI,DJ,DY,DJR,DJI
1045 .EQ.
IF (NP-1) CALL RSP1(X,NG,NGAUSS,A,EPS,NP,R,DR)
1046 .GE.
IF (NP0) CALL RSP2(X,NG,A,EPS,NP,R,DR)
1047 .EQ.
IF (NP-2) CALL RSP3(X,NG,NGAUSS,A,EPS,R,DR)
1052 V=1Q0/(MRR*MRR+MRI*MRI)
1070 .GT.
IF (NMAXNPN1) PRINT 9000,NMAX,NPN1
1071 .GT.
IF (NMAXNPN1) STOP
1072 9000 FORMAT(' nmax =
',I2,', i.e., greater than
',I3)
1073 TB=TA*QSQRT(MRR*MRR+MRI*MRI)
1074 TB=QMAX1(TB,QFLOAT(NMAX))
1075 NNMAX1=8.0Q0*QSQRT(QMAX1(TA,QFLOAT(NMAX)))+3Q0
1076 NNMAX2=(TB+4Q0*(TB**0.33333Q0)+8.0Q0*QSQRT(TB))
1077 NNMAX2=NNMAX2-NMAX+5
1078 CALL BESS(Z,ZR,ZI,NG,NMAX,NNMAX1,NNMAX2)
1082 C**********************************************************************
1084 C Calculation of the functions R(I)=r(y(I))**2 and *
1085 C DR(I)=((d/dy)r(y))/r(y) and horizontal semi-axis A *
1086 C for a spheroid specified by the parameters REV (equal-volume- *
1087 C sphere radius) and EPS=A/B (ratio of the semi-axes). *
1088 C Y(I)=arccos(X(I)) *
1092 C**********************************************************************
1094 SUBROUTINE RSP1 (X,NG,NGAUSS,REV,EPS,NP,R,DR)
1095 IMPLICIT REAL*16 (A-H,O-Z)
1096 REAL*16 X(NG),R(NG),DR(NG)
1097 A=REV*EPS**(1Q0/3Q0)
1115 C**********************************************************************
1117 C Calculation of the functions R(I)=r(y(i))**2 and *
1118 C DR(I)=((d/dy)r(y))/r(y) and parameter R0 for a Chebyshev *
1119 C particle specified by the parameters REV (equal-volume-sphere *
1120 C radius), EPS, and N. *
1121 C Y(I)=arccos(X(I)) *
1125 C**********************************************************************
1127 SUBROUTINE RSP2 (X,NG,REV,EPS,N,R,DR)
1128 IMPLICIT REAL*16 (A-H,O-Z)
1129 REAL*16 X(NG),R(NG),DR(NG)
1134 A=1Q0+1.5Q0*EP*(DN4-2Q0)/(DN4-1Q0)
1137 .EQ.
IF (IN) A=A-3Q0*EPS*(1Q0+0.25Q0*EP)/
1138 * (DN-1Q0)-0.25Q0*EP*EPS/(9Q0*DN-1Q0)
1139 R0=REV*A**(-1Q0/3Q0)
1142 RI=R0*(1Q0+EPS*QCOS(XI))
1144 DR(I)=-R0*EPS*DNP*QSIN(XI)/RI
1149 C**********************************************************************
1151 C Calculation of the functions R(I)=R(Y(I))**2 and *
1152 C DR(I)=((d/dy)r(y))/r(y) *
1153 C for a cylinder specified by the parameters REV (equal-volume- *
1154 C sphere radius) and EPS=A/H (ratio of radius to semi-length) *
1155 C Y(I)=arccos(X(I)) *
1159 C**********************************************************************
1161 SUBROUTINE RSP3 (X,NG,NGAUSS,REV,EPS,R,DR)
1162 IMPLICIT REAL*16 (A-H,O-Z)
1163 REAL*16 X(NG),R(NG),DR(NG)
1164 H=REV*( (2Q0/(3Q0*EPS*EPS))**(1Q0/3Q0) )
1169 .GT.
IF (SI/COA/H) GO TO 20
1183 C**********************************************************************
1185 C Calculation of spherical Bessel functions of the first kind *
1186 C J(I,N) = j_n(x) and second kind Y(I,N) = y_n(x) *
1187 C of real-valued argument X(I) and first kind JR(I,N)+I*JI(I,N) = *
1188 C = j_n(z) of complex argument Z(I)=XR(I)+I*XI(I), as well as *
1191 C DJ(I,N) = (1/x)(d/dx)(x*j_n(x)) , *
1192 C DY(I,N) = (1/x)(d/dx)(x*y_n(x)) , *
1193 C DJR(I,N) = Re ((1/z)(d/dz)(z*j_n(x)) , *
1194 C DJI(I,N) = Im ((1/z)(d/dz)(z*j_n(x)) . *
1198 C X,XR,XI - arguments *
1200 C Arrays J,Y,JR,JI,DJ,DY,DJR,DJI are in *
1202 C Parameters NNMAX1 and NMAX2 determine computational accuracy *
1204 C**********************************************************************
1206 SUBROUTINE BESS (X,XR,XI,NG,NMAX,NNMAX1,NNMAX2)
1208 IMPLICIT REAL*16 (A-H,O-Z)
1209 REAL*16 X(NG),XR(NG),XI(NG),
1210 * J(NPNG2,NPN1),Y(NPNG2,NPN1),JR(NPNG2,NPN1),
1211 * JI(NPNG2,NPN1),DJ(NPNG2,NPN1),DY(NPNG2,NPN1),
1212 * DJR(NPNG2,NPN1),DJI(NPNG2,NPN1),
1213 * AJ(NPN1),AY(NPN1),AJR(NPN1),AJI(NPN1),
1214 * ADJ(NPN1),ADY(NPN1),ADJR(NPN1),
1216 COMMON /CBESS/ J,Y,JR,JI,DJ,DY,DJR,DJI
1220 CALL RJB(XX,AJ,ADJ,NMAX,NNMAX1)
1221 CALL RYB(XX,AY,ADY,NMAX)
1224 CALL CJB(YR,YI,AJR,AJI,ADJR,ADJI,NMAX,NNMAX2)
1238 C**********************************************************************
1240 C Calculation of spherical Bessel functions of the first kind j *
1241 C of real-valued argument x of orders from 1 to NMAX by using *
1242 C backward recursion. Parametr NNMAX determines numerical accuracy. *
1243 C U - function (1/x)(d/dx)(x*j(x)) *
1245 C**********************************************************************
1247 SUBROUTINE RJB(X,Y,U,NMAX,NNMAX)
1248 IMPLICIT REAL*16 (A-H,O-Z)
1249 REAL*16 Y(NMAX),U(NMAX),Z(900)
1252 Z(L)=1Q0/(QFLOAT(2*L+1)*XX)
1256 Z(I1)=1Q0/(QFLOAT(2*I1+1)*XX-Z(I1+1))
1266 U(I)=YI1-QFLOAT(I)*YI*XX
1272 C**********************************************************************
1274 C Calculation of spherical Bessel functions of the second kind y *
1275 C of real-valued argument x of orders from 1 to NMAX by using upward*
1276 C recursion. V - function (1/x)(d/dx)(x*y(x)) *
1278 C**********************************************************************
1280 SUBROUTINE RYB(X,Y,V,NMAX)
1281 IMPLICIT REAL*16 (A-H,O-Z)
1282 REAL*16 Y(NMAX),V(NMAX)
1290 Y(2)=(-3Q0*X3+X1)*C-3Q0*X2*S
1293 5 Y(I+1)=QFLOAT(2*I+1)*X1*Y(I)-Y(I-1)
1296 10 V(I)=Y(I-1)-QFLOAT(I)*X1*Y(I)
1300 C**********************************************************************
1302 C Calculation of spherical Bessel functions of the first kind *
1303 C j=JR+I*JI of complex argument x=XR+I*XI of orders from 1 to NMAX *
1304 C by using backward recursion. Parametr NNMAX determines numerical *
1305 C accuracy. U=UR+I*UI - function (1/x)(d/dx)(x*j(x)) *
1307 C**********************************************************************
1309 SUBROUTINE CJB (XR,XI,YR,YI,UR,UI,NMAX,NNMAX)
1311 IMPLICIT REAL*16 (A-H,O-Z)
1312 REAL*16 YR(NMAX),YI(NMAX),UR(NMAX),UI(NMAX)
1313 REAL*16 CYR(NPN1),CYI(NPN1),CZR(2500),CZI(2500),
1314 * CUR(NPN1),CUI(NPN1)
1316 XRXI=1Q0/(XR*XR+XI*XI)
1319 QF=1Q0/QFLOAT(2*L+1)
1326 AR=QF*CXXR-CZR(I1+1)
1327 AI=QF*CXXI-CZI(I1+1)
1328 ARI=1Q0/(AR*AR+AI*AI)
1334 ARI=1Q0/(AR*AR+AI*AI)
1337 CR=QCOS(XR)*QCOSH(XI)
1338 CI=-QSIN(XR)*QSINH(XI)
1341 CY0R=AR*CXXR-AI*CXXI
1342 CY0I=AI*CXXR+AR*CXXI
1343 CY1R=CY0R*CZR(1)-CY0I*CZI(1)
1344 CY1I=CY0I*CZR(1)+CY0R*CZI(1)
1345 AR=CY1R*CXXR-CY1I*CXXI
1346 AI=CY1I*CXXR+CY1R*CXXI
1361 CYIR=CYI1R*CZR(I)-CYI1I*CZI(I)
1362 CYII=CYI1I*CZR(I)+CYI1R*CZI(I)
1363 AR=CYIR*CXXR-CYII*CXXI
1364 AI=CYII*CXXR+CYIR*CXXI
1379 C**********************************************************************
1381 C calculation of the T(0) matrix for an axially symmetric particle *
1383 C Output information: *
1385 C Arrays TR1 and TI1 (real and imaginary parts of the *
1386 C T(0) matrix) are in COMMON /CT/ *
1388 C**********************************************************************
1390 SUBROUTINE TMATR0 (NGAUSS,X,W,AN,ANN,S,SS,PPI,PIR,PII,R,DR,DDR,
1391 * DRR,DRI,NMAX,NCHECK)
1393 IMPLICIT REAL*16 (A-H,O-Z)
1394 REAL*16 X(NPNG2),W(NPNG2),AN(NPN1),S(NPNG2),SS(NPNG2),
1395 * R(NPNG2),DR(NPNG2),SIG(NPN2),
1396 * J(NPNG2,NPN1),Y(NPNG2,NPN1),
1397 * JR(NPNG2,NPN1),JI(NPNG2,NPN1),DJ(NPNG2,NPN1),
1398 * DY(NPNG2,NPN1),DJR(NPNG2,NPN1),
1399 * DJI(NPNG2,NPN1),DDR(NPNG2),DRR(NPNG2),
1400 * D1(NPNG2,NPN1),D2(NPNG2,NPN1),
1401 * DRI(NPNG2),DS(NPNG2),DSS(NPNG2),RR(NPNG2),
1402 * DV1(NPN1),DV2(NPN1)
1404 REAL*16 R11(NPN1,NPN1),R12(NPN1,NPN1),
1405 * R21(NPN1,NPN1),R22(NPN1,NPN1),
1406 * I11(NPN1,NPN1),I12(NPN1,NPN1),
1407 * I21(NPN1,NPN1),I22(NPN1,NPN1),
1408 * RG11(NPN1,NPN1),RG12(NPN1,NPN1),
1409 * RG21(NPN1,NPN1),RG22(NPN1,NPN1),
1410 * IG11(NPN1,NPN1),IG12(NPN1,NPN1),
1411 * IG21(NPN1,NPN1),IG22(NPN1,NPN1),
1413 * QR(NPN2,NPN2),QI(NPN2,NPN2),
1414 * RGQR(NPN2,NPN2),RGQI(NPN2,NPN2),
1415 * TQR(NPN2,NPN2),TQI(NPN2,NPN2),
1416 * TRGQR(NPN2,NPN2),TRGQI(NPN2,NPN2)
1417 REAL*8 TR1(NPN2,NPN2),TI1(NPN2,NPN2)
1418 REAL*4 PLUS(NPN6*NPN4*NPN4*8)
1420 & R11,R12,R21,R22,I11,I12,I21,I22,RG11,RG12,RG21,RG22,
1421 & IG11,IG12,IG21,IG22
1422 COMMON /CBESS/ J,Y,JR,JI,DJ,DY,DJR,DJI
1424 COMMON /CTT/ QR,QI,RGQR,RGQI
1430 .EQ.
IF (NCHECK1) THEN
1444 CALL VIG (X(I1),NMAX,0,DV1,DV2)
1470 .EQ..AND..LT.
IF (NCHECK1SIG(N1+N2)0Q0) GO TO 205
1513 C5R=C1R*DRRI-C1I*DRII
1514 C5I=C1I*DRRI+C1R*DRII
1515 B5R=B1R*DRRI-B1I*DRII
1516 B5I=B1I*DRRI+B1R*DRII
1523 AR12=AR12+F1*B2R+F2*B3R
1524 AI12=AI12+F1*B2I+F2*B3I
1525 GR12=GR12+F1*C2R+F2*C3R
1526 GI12=GI12+F1*C2I+F2*C3I
1529 AR21=AR21+F1*B4R+F2*B5R
1530 AI21=AI21+F1*B4I+F2*B5I
1531 GR21=GR21+F1*C4R+F2*C5R
1532 GI21=GI21+F1*C4I+F2*C5I
1535 205 AN12=ANN(N1,N2)*FACTOR
1536 R12(N1,N2)=AR12*AN12
1537 R21(N1,N2)=AR21*AN12
1538 I12(N1,N2)=AI12*AN12
1539 I21(N1,N2)=AI21*AN12
1540 RG12(N1,N2)=GR12*AN12
1541 RG21(N1,N2)=GR21*AN12
1542 IG12(N1,N2)=GI12*AN12
1543 IG21(N1,N2)=GI21*AN12
1568 TQR(K1,K2)=TPIR*TAR21-TPII*TAI21+TPPI*TAR12
1569 TQI(K1,K2)=TPIR*TAI21+TPII*TAR21+TPPI*TAI12
1570 TRGQR(K1,K2)=TPIR*TGR21-TPII*TGI21+TPPI*TGR12
1571 TRGQI(K1,K2)=TPIR*TGI21+TPII*TGR21+TPPI*TGI12
1583 TQR(KK1,KK2)=TPIR*TAR12-TPII*TAI12+TPPI*TAR21
1584 TQI(KK1,KK2)=TPIR*TAI12+TPII*TAR12+TPPI*TAI21
1585 TRGQR(KK1,KK2)=TPIR*TGR12-TPII*TGI12+TPPI*TGR21
1586 TRGQI(KK1,KK2)=TPIR*TGI12+TPII*TGR12+TPPI*TGI21
1592 QR(N1,N2)=TQR(N1,N2)
1593 QI(N1,N2)=TQI(N1,N2)
1594 RGQR(N1,N2)=TRGQR(N1,N2)
1595 RGQI(N1,N2)=TRGQI(N1,N2)
1597 CALL TT(NMAX,NCHECK)
1601 C**********************************************************************
1603 .GE.
C Calculation of the T(M) matrix, M1, for an axially symmetric *
1606 C Input parameters: *
1609 C NG = NGAUSS*2 - number of gaussian division points on the *
1611 C W - quadrature weights *
1612 C AN,ANN - see subroutine CONST *
1613 C S,SS - see subroutine CONST *
1614 C ARRAYS DV1,DV2,DV3,DV4 are in COMMON /DV/ - *
1615 C see subroutine DVIG *
1616 C PPI,PIR,PII - see subroutine VARY *
1617 C R J DR - see subroutines RSP1 and RSP2 *
1618 C DDR,DRR,DRI - see subroutine VARY *
1619 C NMAX - dimension of the T(M) matrix *
1620 C Arrays J,Y,JR,JI,DJ,DY,DJR,DJI are in *
1621 C COMMON /CBESS/ - see subroutine BESS *
1623 C Output parameters: *
1625 C Arrays TR1,TI1 (real and imaginary parts of the T(M) matrix) *
1626 C are in COMMON /CT/ *
1628 C**********************************************************************
1630 SUBROUTINE TMATR (M,NGAUSS,X,W,AN,ANN,S,SS,PPI,PIR,PII,R,DR,DDR,
1631 * DRR,DRI,NMAX,NCHECK)
1633 IMPLICIT REAL*16 (A-H,O-Z)
1634 REAL*16 X(NPNG2),W(NPNG2),AN(NPN1),S(NPNG2),SS(NPNG2),
1635 * R(NPNG2),DR(NPNG2),SIG(NPN2),
1636 * J(NPNG2,NPN1),Y(NPNG2,NPN1),
1637 * JR(NPNG2,NPN1),JI(NPNG2,NPN1),DJ(NPNG2,NPN1),
1638 * DY(NPNG2,NPN1),DJR(NPNG2,NPN1),
1639 * DJI(NPNG2,NPN1),DDR(NPNG2),DRR(NPNG2),
1640 * D1(NPNG2,NPN1),D2(NPNG2,NPN1),
1641 * DRI(NPNG2),DS(NPNG2),DSS(NPNG2),RR(NPNG2),
1642 * DV1(NPN1),DV2(NPN1)
1643 REAL*16 R11(NPN1,NPN1),R12(NPN1,NPN1),
1644 * R21(NPN1,NPN1),R22(NPN1,NPN1),
1645 * I11(NPN1,NPN1),I12(NPN1,NPN1),
1646 * I21(NPN1,NPN1),I22(NPN1,NPN1),
1647 * RG11(NPN1,NPN1),RG12(NPN1,NPN1),
1648 * RG21(NPN1,NPN1),RG22(NPN1,NPN1),
1649 * IG11(NPN1,NPN1),IG12(NPN1,NPN1),
1650 * IG21(NPN1,NPN1),IG22(NPN1,NPN1),
1652 * QR(NPN2,NPN2),QI(NPN2,NPN2),
1653 * RGQR(NPN2,NPN2),RGQI(NPN2,NPN2),
1654 * TQR(NPN2,NPN2),TQI(NPN2,NPN2),
1655 * TRGQR(NPN2,NPN2),TRGQI(NPN2,NPN2)
1656 REAL*8 TR1(NPN2,NPN2),TI1(NPN2,NPN2)
1657 REAL*4 PLUS(NPN6*NPN4*NPN4*8)
1659 & R11,R12,R21,R22,I11,I12,I21,I22,RG11,RG12,RG21,RG22,
1660 & IG11,IG12,IG21,IG22
1661 COMMON /CBESS/ J,Y,JR,JI,DJ,DY,DJR,DJI
1663 COMMON /CTT/ QR,QI,RGQR,RGQI
1670 .EQ.
IF (NCHECK1) THEN
1685 CALL VIG (X(I1),NMAX,M,DV1,DV2)
1767 C5R=C1R*DRRI-C1I*DRII
1768 C5I=C1I*DRRI+C1R*DRII
1769 B5R=B1R*DRRI-B1I*DRII
1770 B5I=B1I*DRRI+B1R*DRII
1782 C8R=C2R*DRRI-C2I*DRII
1783 C8I=C2I*DRRI+C2R*DRII
1784 B8R=B2R*DRRI-B2I*DRII
1785 B8I=B2I*DRRI+B2R*DRII
1792 .EQ..AND..GT.
IF (NCHECK1SI0Q0) GO TO 150
1799 .EQ.
IF (NCHECK1) GO TO 160
1803 AR12=AR12+F1*B2R+F2*B3R
1804 AI12=AI12+F1*B2I+F2*B3I
1805 GR12=GR12+F1*C2R+F2*C3R
1806 GI12=GI12+F1*C2I+F2*C3I
1809 AR21=AR21+F1*B4R+F2*B5R
1810 AI21=AI21+F1*B4I+F2*B5I
1811 GR21=GR21+F1*C4R+F2*C5R
1812 GI21=GI21+F1*C4I+F2*C5I
1813 .EQ.
IF (NCHECK1) GO TO 200
1818 AR22=AR22+E1*B6R+E2*B7R+E3*B8R
1819 AI22=AI22+E1*B6I+E2*B7I+E3*B8I
1820 GR22=GR22+E1*C6R+E2*C7R+E3*C8R
1821 GI22=GI22+E1*C6I+E2*C7I+E3*C8I
1823 AN12=ANN(N1,N2)*FACTOR
1824 R11(N1,N2)=AR11*AN12
1825 R12(N1,N2)=AR12*AN12
1826 R21(N1,N2)=AR21*AN12
1827 R22(N1,N2)=AR22*AN12
1828 I11(N1,N2)=AI11*AN12
1829 I12(N1,N2)=AI12*AN12
1830 I21(N1,N2)=AI21*AN12
1831 I22(N1,N2)=AI22*AN12
1832 RG11(N1,N2)=GR11*AN12
1833 RG12(N1,N2)=GR12*AN12
1834 RG21(N1,N2)=GR21*AN12
1835 RG22(N1,N2)=GR22*AN12
1836 IG11(N1,N2)=GI11*AN12
1837 IG12(N1,N2)=GI12*AN12
1838 IG21(N1,N2)=GI21*AN12
1839 IG22(N1,N2)=GI22*AN12
1873 TQR(K1,K2)=TPIR*TAR21-TPII*TAI21+TPPI*TAR12
1874 TQI(K1,K2)=TPIR*TAI21+TPII*TAR21+TPPI*TAI12
1875 TRGQR(K1,K2)=TPIR*TGR21-TPII*TGI21+TPPI*TGR12
1876 TRGQI(K1,K2)=TPIR*TGI21+TPII*TGR21+TPPI*TGI12
1878 TQR(K1,KK2)=TPIR*TAR11-TPII*TAI11+TPPI*TAR22
1879 TQI(K1,KK2)=TPIR*TAI11+TPII*TAR11+TPPI*TAI22
1880 TRGQR(K1,KK2)=TPIR*TGR11-TPII*TGI11+TPPI*TGR22
1881 TRGQI(K1,KK2)=TPIR*TGI11+TPII*TGR11+TPPI*TGI22
1883 TQR(KK1,K2)=TPIR*TAR22-TPII*TAI22+TPPI*TAR11
1884 TQI(KK1,K2)=TPIR*TAI22+TPII*TAR22+TPPI*TAI11
1885 TRGQR(KK1,K2)=TPIR*TGR22-TPII*TGI22+TPPI*TGR11
1886 TRGQI(KK1,K2)=TPIR*TGI22+TPII*TGR22+TPPI*TGI11
1888 TQR(KK1,KK2)=TPIR*TAR12-TPII*TAI12+TPPI*TAR21
1889 TQI(KK1,KK2)=TPIR*TAI12+TPII*TAR12+TPPI*TAI21
1890 TRGQR(KK1,KK2)=TPIR*TGR12-TPII*TGI12+TPPI*TGR21
1891 TRGQI(KK1,KK2)=TPIR*TGI12+TPII*TGR12+TPPI*TGI21
1897 QR(N1,N2)=TQR(N1,N2)
1898 QI(N1,N2)=TQI(N1,N2)
1899 RGQR(N1,N2)=TRGQR(N1,N2)
1900 RGQI(N1,N2)=TRGQI(N1,N2)
1906 C*****************************************************************
1908 c Calculation of the functiONS
1909 c DV1(n)=dvig(0,m,n,arccos x)
1911 c DV2(n)=[d/d(arccos x)] dvig(0,m,n,arccos x)
1915 SUBROUTINE VIG (X,NMAX,M,DV1,DV2)
1917 IMPLICIT REAL*16 (A-H,O-Z)
1918 REAL*16 DV1(NPN1), DV2(NPN1)
1926 .NE.
IF (M0) GO TO 20
1933 D3=(QN2*X*D2-QN*D1)/QN1
1934 DER=QS1*(QN1*QN/QN2)*(-D1+D3)
1944 A=A*QSQRT(QFLOAT(I2-1)/QFLOAT(I2))*QS
1952 QNM=QSQRT(QN*QN-QMM)
1953 QNM1=QSQRT(QN1*QN1-QMM)
1954 D3=(QN2*X*D2-QNM*D1)/QNM1
1955 DER=QS1*(-QN1*QNM*D1+QN*QNM1*D3)/QN2
1964 C**********************************************************************
1966 C Calculation of the matrix T = - RG(Q) * (Q**(-1)) *
1968 C Input infortmation is in COMMON /CTT/ *
1969 C Output information is in COMMON /CT/ *
1971 C**********************************************************************
1973 SUBROUTINE TT(NMAX,NCHECK)
1975 IMPLICIT REAL*8 (A-H,O-Z)
1976 REAL*16 F(NPN2,NPN2),B(NPN2),WORK(NPN2),COND,
1977 * QR(NPN2,NPN2),QI(NPN2,NPN2),
1978 * RGQR(NPN2,NPN2),RGQI(NPN2,NPN2),
1979 * A(NPN2,NPN2),C(NPN2,NPN2),D(NPN2,NPN2),E(NPN2,NPN2)
1980 REAL*8 TR1(NPN2,NPN2),TI1(NPN2,NPN2)
1981 COMPLEX*32 ZQ(NPN2,NPN2),ZW(NPN2)
1982 INTEGER IPIV(NPN2),IPVT(NPN2)
1984 COMMON /CTT/ QR,QI,RGQR,RGQI
1988 C Inversion from LAPACK
1992 ZQ(I,J)=QCMPLX(QR(I,J),QI(I,J))
1996 CALL ZGETRF(NNMAX,NNMAX,ZQ,NPN2,IPIV,INFO)
1997 .NE.
IF (INFO0) WRITE (6,1100) INFO
1998 CALL ZGETRI(NNMAX,ZQ,NPN2,IPIV,ZW,NPN2,INFO)
1999 .NE.
IF (INFO0) WRITE (6,1100) INFO
2001 1100 FORMAT ('warning: info=
', I2)
2021 C********************************************************************
2023 C Calculation of the expansion coefficients for the (I,Q,U,V) *
2026 C Input parameters: *
2028 C LAM - wavelength of light *
2029 C CSCA - scattering cross section *
2030 C TR and TI - elements of the t-matrix. Transferred through *
2032 C NMAX - dimension of T(M) matrices *
2034 C Output infortmation: *
2036 C ALF1,...,ALF4,BET1,BET2 - expansion coefficients *
2037 C LMAX - number of coefficients minus 1 *
2039 C********************************************************************
2041 SUBROUTINE GSP(NMAX,CSCA,LAM,ALF1,ALF2,ALF3,ALF4,BET1,BET2,LMAX)
2043 IMPLICIT REAL*8 (A-B,D-H,O-Z),COMPLEX*16 (C)
2046 REAL*8 CSCA,SSI(NPL),SSJ(NPN1),
2047 & ALF1(NPL),ALF2(NPL),ALF3(NPL),
2048 & ALF4(NPL),BET1(NPL),BET2(NPL),
2049 & TR1(NPL1,NPN4),TR2(NPL1,NPN4),
2050 & TI1(NPL1,NPN4),TI2(NPL1,NPN4),
2051 & G1(NPL1,NPN6),G2(NPL1,NPN6),
2052 & AR1(NPN4),AR2(NPN4),AI1(NPN4),AI2(NPN4),
2053 & FR(NPN4,NPN4),FI(NPN4,NPN4),FF(NPN4,NPN4)
2054 REAL*4 B1R(NPL1,NPL1,NPN4),B1I(NPL1,NPL1,NPN4),
2055 & B2R(NPL1,NPL1,NPN4),B2I(NPL1,NPL1,NPN4),
2056 & D1(NPL1,NPN4,NPN4),D2(NPL1,NPN4,NPN4),
2057 & D3(NPL1,NPN4,NPN4),D4(NPL1,NPN4,NPN4),
2058 & D5R(NPL1,NPN4,NPN4),D5I(NPL1,NPN4,NPN4),
2059 & PLUS1(NPN6*NPN4*NPN4*8)
2061 & TR11(NPN6,NPN4,NPN4),TR12(NPN6,NPN4,NPN4),
2062 & TR21(NPN6,NPN4,NPN4),TR22(NPN6,NPN4,NPN4),
2063 & TI11(NPN6,NPN4,NPN4),TI12(NPN6,NPN4,NPN4),
2064 & TI21(NPN6,NPN4,NPN4),TI22(NPN6,NPN4,NPN4)
2065 COMPLEX*16 CIM(NPN1)
2067 COMMON /TMAT/ TR11,TR12,TR21,TR22,TI11,TI12,TI21,TI22
2068 COMMON /CBESS/ B1R,B1I,B2R,B2I
2070 EQUIVALENCE ( PLUS1(1),TR11(1,1,1) )
2071 EQUIVALENCE (D1(1,1,1),PLUS1(1)),
2072 & (D2(1,1,1),PLUS1(NPL1*NPN4*NPN4+1)),
2073 & (D3(1,1,1),PLUS1(NPL1*NPN4*NPN4*2+1)),
2074 & (D4(1,1,1),PLUS1(NPL1*NPN4*NPN4*3+1)),
2075 & (D5R(1,1,1),PLUS1(NPL1*NPN4*NPN4*4+1))
2091 .LE.
IF(INMAX) SSJ(I)=DSQRT(SI)
2107 C ***** CALCULATION OF THE ARRAYS B1 AND B2 *****
2117 3300 FORMAT (' b1
and b2
')
2120 C ***** CALCULATION OF THE ARRAYS T1 AND T2 *****
2157 C ***** END OF THE CALCULATION OF THE ARRAYS T1 AND T2 *****
2163 C ***** CALCULATION OF THE ARRAYS A1 AND A2 *****
2165 CALL CCG(N,N1,NMAX,K1,K2,G1)
2166 NNMAX=MIN0(NMAX,N1+N)
2167 NNMIN=MAX0(1,IABS(N-N1))
2169 DO 15 NN=NNMIN,NNMAX
2172 M1MAX=MIN0(N,NN)+NPN6
2184 .EQ.
IF(M0) GO TO 12
2186 RR1=RR1+TR1(M2,NN)*SIG
2187 RI1=RI1+TI1(M2,NN)*SIG
2188 RR2=RR2+TR2(M2,NN)*SIG
2189 RI2=RI2+TI2(M2,NN)*SIG
2190 12 AAR1=AAR1+SSS*RR1
2197 AR1(NN)=AAR1*XR-AAI1*XI
2198 AI1(NN)=AAR1*XI+AAI1*XR
2199 AR2(NN)=AAR2*XR-AAI2*XI
2200 AI2(NN)=AAR2*XI+AAI2*XR
2203 C ***** END OF THE CALCULATION OF THE ARRAYS A1 AND A2 ****
2205 CALL CCG(N,N1,NMAX,K3,K4,G2)
2210 DO 30 M1=M1MIN,M1MAX
2215 DO 25 NN=NNMIN,NNMAX
2218 BBR1=BBR1+SSS*AR1(NN)
2219 BBI1=BBI1+SSS*AI1(NN)
2220 BBR2=BBR2+SSS*AR2(NN)
2221 BBI2=BBI2+SSS*AI2(NN)
2231 C ***** END OF THE CALCULATION OF THE ARRAYS B1 AND B2 ****
2233 C ***** CALCULATION OF THE ARRAYS D1,D2,D3,D4, AND D5 *****
2236 3301 FORMAT(' d1, d2, ...
')
2242 NN1MAX=NMAX1+MIN0(N,NN)
2243 DO 180 M1=M1MIN,M1MAX
2248 DO 150 NN1=NN1MIN,NN1MAX
2258 DD1=DD1+XX*(X1*X3+X2*X4)
2259 DD2=DD2+XX*(X5*X7+X6*X8)
2268 DO 186 M1=M1MIN,M1MAX
2276 DO 183 NN1=NN1MIN,NN1MAX
2286 DD3=DD3+XX*(X1*X5+X2*X6)
2287 DD4=DD4+XX*(X3*X7+X4*X8)
2288 DD5R=DD5R+XX*(X3*X5+X4*X6)
2289 DD5I=DD5I+XX*(X4*X5-X3*X6)
2299 C ***** END OF THE CALCULATION OF THE D-ARRAYS *****
2301 C ***** CALCULATION OF THE EXPANSION COEFFICIENTS *****
2304 3303 FORMAT (' g1, g2, ...
')
2306 DK=LAM*LAM/(4D0*CSCA*DACOS(-1D0))
2317 NNMIN=MAX0(1,IABS(N-L))
2318 NNMAX=MIN0(NMAX,N+L)
2319 .LT.
IF(NNMAXNNMIN) GO TO 290
2320 CALL CCG(N,L,NMAX,K1,K2,G1)
2321 .GE.
IF(L2) CALL CCG(N,L,NMAX,K5,K6,G2)
2323 DO 280 NN=NNMIN,NNMAX
2331 DO 270 M1=M1MIN,M1MAX
2333 .GE.
IF(M0) SSS1=G1(M1,NNN)
2334 .LT.
IF(M0) SSS1=G1(NPN6-M,NNN)*SI
2335 DM1=DM1+SSS1*D1(M1,NN,N)
2336 DM2=DM2+SSS1*D2(M1,NN,N)
2339 SSS=G1(NPN6+1,NNN)*FFN
2342 .LT.
IF(L2) GO TO 280
2351 DO 275 M1=M1MIN,M1MAX
2354 DM3=DM3+SSS1*D3(M1,NN,N)
2355 DM4=DM4+SSS1*D4(M1,NN,N)
2356 DM5R=DM5R+SSS1*D5R(M1,NN,N)
2357 DM5I=DM5I+SSS1*D5I(M1,NN,N)
2361 SSS=G2(NPN4,NNN)*FFN
2379 .LT.
IF(DABS(G1L)1D-6) GO TO 500
2385 C****************************************************************
2387 C Calculation of the quantities F(N+1)=0.5*ln(n!)
2397 F(I)=F(I1)+0.5D0*DLOG(DFLOAT(I1))
2402 C************************************************************
2404 C Calculation of the array SSIGN(N+1)=sign(n)
2412 SSIGN(N)=-SSIGN(N-1)
2417 C******************************************************************
2419 C Calculation of Clebsch-Gordan coefficients
2421 C for given n and n1. m1=mm-m, index mm is found from m as
2424 C Input parameters : N,N1,NMAX,K1,K2
2429 C Output parameters : GG(M+NPN6,NN+1) - array of the corresponding
2432 .LE.
C /M1/=/M*(K1-1)+K2/N1
2433 .LE.
C NNMIN(N+N1,NMAX)
2434 .GE.
C NNMAX(/MM/,/N-N1/)
2435 .LE..LE.
C If K1=1 and K2=0, then 0MN
2437 SUBROUTINE CCG(N,N1,NMAX,K1,K2,GG)
2439 IMPLICIT REAL*8 (A-H,O-Z)
2440 REAL*8 GG(NPL1,NPN6),CD(0:NPN5),CU(0:NPN5)
2443 .LE.
& AND.N1NMAX+N.
2445 .LE.
& AND.NNMAX) GO TO 1
2448 5001 FORMAT(' error in
SUBROUTINE ccg')
2449 1 NNF=MIN0(N+N1,NMAX)
2452 .EQ..AND..EQ.
IF(K11K20) MIN=NPN6
2457 .GT.
IF(IABS(M1)N1) GO TO 90
2458 NNL=MAX0(IABS(MM),IABS(N-N1))
2459 .GT.
IF(NNLNNF) GO TO 90
2462 .EQ.
IF (NNLNNU) NNM=NNL
2463 CALL CCGIN(N,N1,M,MM,C)
2465 .EQ.
IF (NNLNNF) GO TO 50
2468 DO 7 NN=NNL+1,MIN0(NNM,NNF)
2469 A=DFLOAT((NN+MM)*(NN-MM)*(N1-N+NN))
2470 A=A*DFLOAT((N-N1+NN)*(N+N1-NN+1)*(N+N1+NN+1))
2472 A=A*DFLOAT((2*NN+1)*(2*NN-1))
2474 B=0.5D0*DFLOAT(M-M1)
2476 .EQ.
IF(NN1) GO TO 5
2477 B=DFLOAT(2*NN*(NN-1))
2478 B=DFLOAT((2*M-MM)*NN*(NN-1)-MM*N*(N+1)+
2480 D=DFLOAT(4*(NN-1)*(NN-1))
2481 D=D*DFLOAT((2*NN-3)*(2*NN-1))
2482 D=DFLOAT((NN-MM-1)*(NN+MM-1)*(N1-N+NN-1))/D
2483 D=D*DFLOAT((N-N1+NN-1)*(N+N1-NN+2)*(N+N1+NN))
2490 .LE.
IF (NNFNNM) GO TO 50
2491 CALL DIRECT(N,M,N1,M1,NNU,MM,C)
2493 .EQ.
IF (NNUNNM+1) GO TO 50
2496 DO 12 NN=NNU-1,NNM+1,-1
2497 A=DFLOAT((NN-MM+1)*(NN+MM+1)*(N1-N+NN+1))
2498 A=A*DFLOAT((N-N1+NN+1)*(N+N1-NN)*(N+N1+NN+2))
2499 A=DFLOAT(4*(NN+1)*(NN+1))/A
2500 A=A*DFLOAT((2*NN+1)*(2*NN+3))
2502 B=DFLOAT(2*(NN+2)*(NN+1))
2503 B=DFLOAT((2*M-MM)*(NN+2)*(NN+1)-MM*N*(N+1)
2505 D=DFLOAT(4*(NN+2)*(NN+2))
2506 D=D*DFLOAT((2*NN+5)*(2*NN+3))
2507 D=DFLOAT((NN+MM+2)*(NN-MM+2)*(N1-N+NN+2))/D
2508 D=D*DFLOAT((N-N1+NN+2)*(N+N1-NN-1)*(N+N1+NN+3))
2516 .LE.
IF (NNNNM) GG(MIND,NN+1)=CU(NN)
2517 .GT.
IF (NNNNM) GG(MIND,NN+1)=CD(NN)
2518 c WRITE (6,*) N,M,N1,M1,NN,MM,GG(MIND,NN+1)
2525 C*********************************************************************
2527 SUBROUTINE DIRECT (N,M,N1,M1,NN,MM,C)
2528 IMPLICIT REAL*8 (A-H,O-Z)
2531 C=F(2*N+1)+F(2*N1+1)+F(N+N1+M+M1+1)+F(N+N1-M-M1+1)
2532 C=C-F(2*(N+N1)+1)-F(N+M+1)-F(N-M+1)-F(N1+M1+1)-F(N1-M1+1)
2537 C*********************************************************************
2539 C Calculation of the Clebcsh-Gordan coefficients
2540 C G=(n,m:n1,mm-m/nn,mm)
2541 C for given n,n1,m,mm, where NN=MAX(/MM/,/N-N1/)
2546 SUBROUTINE CCGIN(N,N1,M,MM,G)
2547 IMPLICIT REAL*8 (A-H,O-Z)
2548 REAL*8 F(900),SSIGN(900)
2553 .GE.
& AND.N1IABS(M1).
2554 .LE.
& AND.IABS(MM)(N+N1)) GO TO 1
2557 5001 FORMAT(' error in subroutine
ccgin')
2558 .GT.
1 IF (IABS(MM)IABS(N-N1)) GO TO 100
2562 .LE.
IF(N1N) GO TO 50
2574 & *DEXP(F(N+M+1)+F(N-M+1)+F(N12+1)+F(N2-N12+2)-F(N2+2)
2575 & -F(N1+M1+1)-F(N1-M1+1)-F(N-N1+MM+1)-F(N-N1-MM+1))
2583 .GE.
IF(MM0) GO TO 150
2588 150 G=A*SSIGN(N+M+1)
2589 & *DEXP(F(2*MM+2)+F(N+N1-MM+1)+F(N+M+1)+F(N1+M1+1)
2590 & -F(N+N1+MM+2)-F(N-N1+MM+1)-F(-N+N1+MM+1)-F(N-M+1)
2597 C*****************************************************************
2599 SUBROUTINE SAREA (D,RAT)
2600 IMPLICIT REAL*8 (A-H,O-Z)
2601 .GE.
IF (D1) GO TO 10
2603 R=0.5D0*(D**(2D0/3D0) + D**(-1D0/3D0)*DASIN(E)/E)
2607 10 E=DSQRT(1D0-1D0/(D*D))
2608 R=0.25D0*(2D0*D**(2D0/3D0) + D**(-4D0/3D0)*DLOG((1D0+E)/(1D0-E))
2615 c****************************************************************
2617 SUBROUTINE SURFCH (N,E,RAT)
2618 IMPLICIT REAL*8 (A-H,O-Z)
2624 CALL GAUSS (NG,0,0,X,W)
2637 S=S+W(I)*A*DSQRT(A2+ENS*ENS)
2638 V=V+W(I)*(DS*A+XI*ENS)*DS*A2
2641 RV=(V*3D0/4D0)**(1D0/3D0)
2646 C********************************************************************
2648 SUBROUTINE SAREAC (EPS,RAT)
2649 IMPLICIT REAL*8 (A-H,O-Z)
2650 RAT=(1.5D0/EPS)**(1D0/3D0)
2651 RAT=RAT/DSQRT( (EPS+2D0)/(2D0*EPS) )
2655 C********************************************************************
2657 C Computation of R1 and R2 for a power law size distribution with
2658 C effective radius A and effective variance B
2660 SUBROUTINE POWER (A,B,R1,R2)
2661 IMPLICIT REAL*8 (A-H,O-Z)
2668 R1=ZEROIN (AX,BX,F,0D0)
2673 C***********************************************************************
2675 DOUBLE PRECISION FUNCTION ZEROIN (AX,BX,F,TOL)
2676 IMPLICIT REAL*8 (A-H,O-Z)
2680 .GT.
IF (TOL11D0) GO TO 10
2689 .GE.
30 IF (DABS(FC)DABS(FB)) GO TO 40
2696 40 TOL1=2D0*EPS*DABS(B)+0.5D0*TOL
2698 .LE.
IF (DABS(XM)TOL1) GO TO 90
2699 .EQ.
44 IF (FB0D0) GO TO 90
2700 .LT.
45 IF (DABS(E)TOL1) GO TO 70
2701 .LE.
46 IF (DABS(FA)DABS(FB)) GO TO 70
2702 .NE.
47 IF (AC) GO TO 50
2710 P=S*(2D0*XM*Q*(Q-R)-(B-A)*(R-1D0))
2711 Q=(Q-1D0)*(R-1D0)*(S-1D0)
2712 .GT.
60 IF (P0D0) Q=-Q
2714 .GE.
IF ((2D0*P)(3D0*XM*Q-DABS(TOL1*Q))) GO TO 70
2715 .GE.
64 IF (PDABS(0.5D0*E*Q)) GO TO 70
2723 .GT.
IF (DABS(D)TOL1) B=B+D
2724 .LE.
IF (DABS(D)TOL1) B=B+DSIGN(TOL1,XM)
2726 .GT.
IF ((FB*(FC/DABS(FC)))0D0) GO TO 20
2732 C***********************************************************************
2734 DOUBLE PRECISION FUNCTION F(R1)
2735 IMPLICIT REAL*8 (A-H,O-Z)
2738 F=(R2-R1)/DLOG(R2/R1)-A
2742 C**********************************************************************
2743 C CALCULATION OF POINTS AND WEIGHTS OF GAUSSIAN QUADRATURE *
2744 C FORMULA. IF IND1 = 0 - ON INTERVAL (-1,1), IF IND1 = 1 - ON *
2745 C INTERVAL (0,1). IF IND2 = 1 RESULTS ARE PRINTED. *
2746 C N - NUMBER OF POINTS *
2747 C Z - DIVISION POINTS *
2749 C**********************************************************************
2751 SUBROUTINE GAUSS ( N,IND1,IND2,Z,W )
2752 IMPLICIT REAL*8 (A-H,P-Z)
2762 .EQ.
IF(I1) X=A-B/((F+A)*F)
2763 .EQ.
IF(I2) X=(Z(N)-A)*4D0+Z(N)
2764 .EQ.
IF(I3) X=(Z(N-1)-Z(N))*1.6D0+Z(N-1)
2765 .GT.
IF(I3) X=(Z(M+1)-Z(M+2))*C+Z(M+3)
2766 .EQ..AND..EQ.
IF(IKIND1) X=0D0
2771 .LE.
IF (NITER100) GO TO 15
2779 20 PC=X*PB+(X*PB-PA)*(DJ-A)/DJ
2783 .GT.
IF(DABS(PB)CHECK*DABS(X)) GO TO 10
2786 .EQ.
IF(IND10) W(M)=B*W(M)
2787 .EQ..AND..EQ.
IF(IKIND1) GO TO 100
2791 .NE.
IF(IND21) GO TO 110
2793 1100 FORMAT(' *** points
and weights of gaussian quadrature formula
',
2794 * ' of
',I4,'-th order
')
2797 105 PRINT 1200,I,ZZ,I,W(I)
2798 1200 FORMAT(' ',4X,'x(',I4,') =
',F17.14,5X,'w(',I4,') =
',F17.14)
2802 1300 FORMAT(' gaussian quadrature formula of
',I4,'-th order is used
')
2804 .EQ.
IF(IND10) GO TO 140
2811 C****************************************************************
2813 SUBROUTINE DISTRB (NNK,YY,WY,NDISTR,AA,BB,GAM,R1,R2,REFF,
2815 IMPLICIT REAL*8 (A-H,O-Z)
2816 REAL*8 YY(NNK),WY(NNK)
2817 .EQ.
IF (NDISTR2) GO TO 100
2818 .EQ.
IF (NDISTR3) GO TO 200
2819 .EQ.
IF (NDISTR4) GO TO 300
2820 .EQ.
IF (NDISTR5) GO TO 360
2821 PRINT 1001,AA,BB,GAM
2822 1001 FORMAT('modified gamma distribution, alpha=
',F6.4,' r_c=
',
2823 & F6.4,' gamma=
',F6.4)
2830 Y=Y*DEXP(-A2*(X**GAM))
2834 100 PRINT 1002,AA,BB
2835 1002 FORMAT('log-normal distribution, r_g=
',F8.4,
2836 & ' [ln(sigma_g)]**2=
',F6.4)
2841 Y=DEXP(-Y*Y*0.5D0/BB)/X
2846 1003 FORMAT('power law distribution of hansen & travis 1974
')
2852 300 PRINT 1004,AA,BB
2853 1004 FORMAT ('gamma distribution, a=
',F8.4,' b=
',F6.4)
2858 X=(X**B2)*DEXP(-X*DAB)
2863 1005 FORMAT ('modified
power law distribution, alpha=
',D10.4)
2866 .LE.
IF (XR1) WY(I)=WY(I)
2867 .GT.
IF (XR1) WY(I)=WY(I)*(X/R1)**BB
2886 REFF=REFF+X*X*X*WY(I)
2893 VEFF=VEFF+XI*XI*X*X*WY(I)
2895 VEFF=VEFF/(G*REFF*REFF)
2899 C*************************************************************
2901 SUBROUTINE HOVENR(L1,A1,A2,A3,A4,B1,B2)
2902 IMPLICIT REAL*8 (A-H,O-Z)
2903 REAL*8 A1(L1),A2(L1),A3(L1),A4(L1),B1(L1),B2(L1)
2907 DL=DFLOAT(LL)*2D0+1D0
2915 .GE..AND..GE.
IF(LL1DABS(AA1)DL) KONTR=2
2916 .GE.
IF(DABS(AA2)DL) KONTR=2
2917 .GE.
IF(DABS(AA3)DL) KONTR=2
2918 .GE.
IF(DABS(AA4)DL) KONTR=2
2919 .GE.
IF(DABS(BB1)DDL) KONTR=2
2920 .GE.
IF(DABS(BB2)DDL) KONTR=2
2921 .EQ.
IF(KONTR2) PRINT 3000,LL
2929 .LE.
IF((DL-C*AA1)*(DL-C*AA2)-CC*BB1*BB1-1D-4) KONTR=2
2930 .LE.
IF((DL-C2)*(DL-C3)+C1-1D-4) KONTR=2
2931 .LE.
IF((DL+C2)*(DL-C3)-C1-1D-4) KONTR=2
2932 .LE.
IF((DL-C2)*(DL+C3)-C1-1D-4) KONTR=2
2933 .EQ.
IF(KONTR2) PRINT 4000,LL,C
2936 .EQ.
IF(KONTR1) PRINT 2000
2937 2000 FORMAT('test of van
der mee & hovenier is satisfied
')
2938 3000 FORMAT('test of van
der mee & hovenier is not satisfied, l=
',I3)
2939 4000 FORMAT('test of van
der mee & hovenier is not satisfied, l=
',I3,
2944 C****************************************************************
2946 C CALCULATION OF THE SCATTERING MATRIX FOR GIVEN EXPANSION
2949 C A1,...,B2 - EXPANSION COEFFICIENTS
2950 C LMAX - NUMBER OF COEFFICIENTS MINUS 1
2951 C N - NUMBER OF SCATTERING ANGLES
2952 C THE CORRESPONDING SCATTERING ANGLES ARE GIVEN BY
2953 C 180*(I-1)/(N-1) (DEGREES), WHERE I NUMBERS THE ANGLES
2955 SUBROUTINE MATR(A1,A2,A3,A4,B1,B2,LMAX,NPNA)
2957 IMPLICIT REAL*8 (A-H,O-Z)
2958 REAL*8 A1(NPL),A2(NPL),A3(NPL),A4(NPL),B1(NPL),B2(NPL)
2967 1001 FORMAT(' ',2X,'s
',6X,'alpha1
',6X,'alpha2
',6X,'alpha3
',
2968 & 6X,'alpha4
',7X,'beta1
',7X,'beta2
')
2971 PRINT 1002,L,A1(L1),A2(L1),A3(L1),A4(L1),B1(L1),B2(L1)
2973 1002 FORMAT(' ',I3,6F12.5)
2978 1003 FORMAT(' ',5X,'<
',8X,'f11
',8X,'f22
',8X,'f33
',
2979 & 8X,'f44
',8X,'f12
',8X,'f34
')
2980 D6=DSQRT(6D0)*0.25D0
2996 PP2=0.25D0*(1D0+U)*(1D0+U)
2997 PP3=0.25D0*(1D0-U)*(1D0-U)
3005 .EQ.
IF(LLMAX) GO TO 350
3007 P=(PL1*U*PP1-DL*P1)/DL1
3010 .LT.
350 IF(L2) GO TO 400
3011 F2=F2+(A2(L1)+A3(L1))*PP2
3012 F3=F3+(A2(L1)-A3(L1))*PP3
3015 .EQ.
IF(LLMAX) GO TO 400
3017 PL3=DFLOAT(L1*(L*L-4))
3018 PL4=1D0/DFLOAT(L*(L1*L1-4))
3019 P=(PL1*(PL2-4D0)*PP2-PL3*P2)*PL4
3022 P=(PL1*(PL2+4D0)*PP3-PL3*P3)*PL4
3025 P=(PL1*U*PP4-DSQRT(DFLOAT(L*L-4))*P4)/DSQRT(DFLOAT(L1*L1-4))
3036 PRINT 1004,TB,F11,F22,F33,F44,F12,F34
3039 1004 FORMAT(' ',F6.2,6F11.4)
