ocssw V2020
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 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 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 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, F)
362  IMPLICIT REAL*8 (a-h,o-z)
363  include 'tmd.par.f'
364  real*8 reff,veff,cextin,cscat,walb,asymm
365  real*8 alpha(npl,4),beta(npl,2),f(npna,4,4)
366 Cf2py intent(in) rat
367 Cf2py intent(in) lam
368 Cf2py intent(in) mrr
369 Cf2py intent(in) mri
370 Cf2py intent(in) eps
371 Cf2py intent(in) np
372 Cf2py intent(in) ddelt
373 Cf2py intent(in) ndgs
374 Cf2py intent(in) npnax
375 Cf2py intent(in) axmax
376 Cf2py intent(in) b
377 Cf2py intent(in) gam
378 Cf2py intent(in) ndistr
379 Cf2py intent(in) nkmax
380 Cf2py intent(in) npna
381 Cf2py intent(in) ncoeff
382 Cf2py intent(out) reff
383 Cf2py intent(out) veff
384 Cf2py intent(out) cextin
385 Cf2py intent(out) cscat
386 Cf2py intent(out) walb
387 Cf2py intent(out) asymm
388 Cf2py intent(out) alpha
389 Cf2py intent(out) beta
390 Cf2py intent(out) f
391 
392  real*8 lam,mrr,mri,x(npng2),w(npng2),s(npng2),ss(npng2),
393  * an(npn1),r(npng2),dr(npng2),
394  * ddr(npng2),drr(npng2),dri(npng2),ann(npn1,npn1)
395  real*8 xg(1000),wg(1000),tr1(npn2,npn2),ti1(npn2,npn2),
396  & alph1(npl),alph2(npl),alph3(npl),alph4(npl),bet1(npl),
397  & bet2(npl),xg1(2000),wg1(2000),
398  & al1(npl),al2(npl),al3(npl),al4(npl),be1(npl),be2(npl)
399  real*4
400  & rt11(npn6,npn4,npn4),rt12(npn6,npn4,npn4),
401  & rt21(npn6,npn4,npn4),rt22(npn6,npn4,npn4),
402  & it11(npn6,npn4,npn4),it12(npn6,npn4,npn4),
403  & it21(npn6,npn4,npn4),it22(npn6,npn4,npn4)
404 
405 
406  COMMON /ct/ tr1,ti1
407  COMMON /tmat/ rt11,rt12,rt21,rt22,it11,it12,it21,it22
408  p=dacos(-1d0)
409 
410 C OPEN FILES *******************************************************
411 
412  OPEN (6,file='tmatr.test')
413  OPEN (10,file='tmatr.write')
414 
415 C INPUT DATA ********************************************************
416 
417 C RAT=0.5 D0
418 C NDISTR=3
419 C AXMAX=1D0
420 C NPNAX=2
421 C B=0.1D0
422 C GAM=0.5D0
423 C NKMAX=5
424 C EPS=2D0
425 C NP=-1
426 C LAM=0.5D0
427 C MRR=1.53 d0
428 C MRI=0.008D0
429 C DDELT=0.001D0
430 C NPNA=19
431 C NDGS=2
432 
433  ncheck=0
434  IF (np.EQ.-1.OR.np.EQ.-2) ncheck=1
435  IF (np.GT.0.AND.(-1)**np.EQ.1) ncheck=1
436  WRITE (6,5454) ncheck
437  5454 FORMAT ('NCHECK=',i1)
438  dax=axmax/npnax
439  IF (abs(rat-1d0).GT.1d-8.AND.np.EQ.-1) CALL sarea (eps,rat)
440  if (abs(rat-1d0).GT.1d-8.AND.np.GE.0) CALL surfch(np,eps,rat)
441  IF (abs(rat-1d0).GT.1d-8.AND.np.EQ.-2) CALL sareac (eps,rat)
442 C PRINT 8000, RAT
443  8000 FORMAT ('RAT=',f8.6)
444  IF(np.EQ.-1.AND.eps.GE.1d0) print 7000,eps
445  IF(np.EQ.-1.AND.eps.LT.1d0) print 7001,eps
446  IF(np.GE.0) print 7100,np,eps
447  IF(np.EQ.-2.AND.eps.GE.1d0) print 7150,eps
448  IF(np.EQ.-2.AND.eps.LT.1d0) print 7151,eps
449  print 7400, lam,mrr,mri
450  print 7200, ddelt
451  7000 FORMAT('RANDOMLY ORIENTED OBLATE SPHEROIDS, A/B=',f11.7)
452  7001 FORMAT('RANDOMLY ORIENTED PROLATE SPHEROIDS, A/B=',f11.7)
453  7100 FORMAT('RANDOMLY ORIENTED CHEBYSHEV PARTICLES, T',
454  & i1,'(',f5.2,')')
455  7150 FORMAT('RANDOMLY ORIENTED OBLATE CYLINDERS, D/L=',f11.7)
456  7151 FORMAT('RANDOMLY ORIENTED PROLATE CYLINDERS, D/L=',f11.7)
457  7200 FORMAT ('ACCURACY OF COMPUTATIONS DDELT = ',d8.2)
458  7400 FORMAT('LAM=',f10.6,3x,'MRR=',d10.4,3x,'MRI=',d10.4)
459  ddelt=0.1d0*ddelt
460  DO 600 iax=1,npnax
461  axi=axmax-dax*dfloat(iax-1)
462  r1=0.89031d0*axi
463  r2=1.56538d0*axi
464  nk=int(axi*nkmax/axmax+2)
465  IF (nk.GT.1000) print 8001,nk
466  IF (nk.GT.1000) stop
467  IF (ndistr.EQ.3) CALL power (axi,b,r1,r2)
468  8001 FORMAT ('NK=',i4,' I.E., IS GREATER THAN 1000. ',
469  & 'EXECUTION TERMINATED.')
470  CALL gauss (nk,0,0,xg,wg)
471  z1=(r2-r1)*0.5d0
472  z2=(r1+r2)*0.5d0
473  z3=r1*0.5d0
474  IF (ndistr.EQ.5) GO TO 3
475  DO i=1,nk
476  xg1(i)=z1*xg(i)+z2
477  wg1(i)=wg(i)*z1
478  ENDDO
479  GO TO 4
480  3 DO i=1,nk
481  xg1(i)=z3*xg(i)+z3
482  wg1(i)=wg(i)*z3
483  ENDDO
484  DO i=nk+1,2*nk
485  ii=i-nk
486  xg1(i)=z1*xg(ii)+z2
487  wg1(i)=wg(ii)*z1
488  ENDDO
489  nk=nk*2
490  4 CALL distrb (nk,xg1,wg1,ndistr,axi,b,gam,r1,r2,
491  & reff,veff,p)
492  print 8002,r1,r2
493  8002 FORMAT('R1=',f10.6,' R2=',f10.6)
494  IF (abs(rat-1d0).LE.1d-6) print 8003, reff,veff
495  IF (abs(rat-1d0).GT.1d-6) print 8004, reff,veff
496  8003 FORMAT('EQUAL-VOLUME-SPHERE REFF=',f8.4,' VEFF=',f7.4)
497  8004 FORMAT('EQUAL-SURFACE-AREA-SPHERE REFF=',f8.4,
498  & ' VEFF=',f7.4)
499  print 7250,nk
500  7250 FORMAT('NUMBER OF GAUSSIAN QUADRATURE POINTS ',
501  & 'IN SIZE AVERAGING =',i4)
502  DO i=1,npl
503  alph1(i)=0d0
504  alph2(i)=0d0
505  alph3(i)=0d0
506  alph4(i)=0d0
507  bet1(i)=0d0
508  bet2(i)=0d0
509  alpha(i,1)=0d0
510  alpha(i,2)=0d0
511  alpha(i,3)=0d0
512  alpha(i,4)=0d0
513  beta(i,1)=0d0
514  beta(i,2)=0d0
515  ENDDO
516  cscat=0d0
517  cextin=0d0
518  l1max=0
519  DO 500 ink=1,nk
520  i=nk-ink+1
521  a=rat*xg1(i)
522  xev=2d0*p*a/lam
523  ixxx=xev+4.05d0*xev**0.333333d0
524  inm1=max0(4,ixxx)
525  IF (inm1.GE.npn1) print 7333, npn1
526  IF (inm1.GE.npn1) stop
527  7333 FORMAT('CONVERGENCE IS NOT OBTAINED FOR NPN1=',i3,
528  & '. EXECUTION TERMINATED')
529  qext1=0d0
530  qsca1=0d0
531  DO 50 nma=inm1,npn1
532  nmax=nma
533  mmax=1
534  ngauss=nmax*ndgs
535  IF (ngauss.GT.npng1) print 7340, ngauss
536  IF (ngauss.GT.npng1) stop
537  7340 FORMAT('NGAUSS =',i3,' I.E. IS GREATER THAN NPNG1.',
538  & ' EXECUTION TERMINATED')
539  7334 FORMAT(' NMAX =', i3,' DC2=',d8.2,' DC1=',d8.2)
540  7335 FORMAT(' NMAX1 =', i3,' DC2=',d8.2,
541  & ' DC1=',d8.2)
542  CALL const(ngauss,nmax,mmax,p,x,w,an,ann,s,ss,np,eps)
543  CALL vary(lam,mrr,mri,a,eps,np,ngauss,x,p,ppi,pir,pii,r,
544  & dr,ddr,drr,dri,nmax)
545  CALL tmatr0 (ngauss,x,w,an,ann,s,ss,ppi,pir,pii,r,dr,
546  & ddr,drr,dri,nmax,ncheck)
547  qext=0d0
548  qsca=0d0
549  DO n=1,nmax
550  n1=n+nmax
551  tr1nn=tr1(n,n)
552  ti1nn=ti1(n,n)
553  tr1nn1=tr1(n1,n1)
554  ti1nn1=ti1(n1,n1)
555  dn1=dfloat(2*n+1)
556  qsca=qsca+dn1*(tr1nn*tr1nn+ti1nn*ti1nn
557  & +tr1nn1*tr1nn1+ti1nn1*ti1nn1)
558  qext=qext+(tr1nn+tr1nn1)*dn1
559  ENDDO
560  dsca=abs((qsca1-qsca)/qsca)
561  dext=abs((qext1-qext)/qext)
562 C PRINT 7334, NMAX,DSCA,DEXT
563  qext1=qext
564  qsca1=qsca
565  nmin=dfloat(nmax)/2d0+1d0
566  DO 10 n=nmin,nmax
567  n1=n+nmax
568  tr1nn=tr1(n,n)
569  ti1nn=ti1(n,n)
570  tr1nn1=tr1(n1,n1)
571  ti1nn1=ti1(n1,n1)
572  dn1=dfloat(2*n+1)
573  dqsca=dn1*(tr1nn*tr1nn+ti1nn*ti1nn
574  & +tr1nn1*tr1nn1+ti1nn1*ti1nn1)
575  dqext=(tr1nn+tr1nn1)*dn1
576  dqsca=abs(dqsca/qsca)
577  dqext=abs(dqext/qext)
578  nmax1=n
579  IF (dqsca.LE.ddelt.AND.dqext.LE.ddelt) GO TO 12
580  10 CONTINUE
581  12 CONTINUE
582 c PRINT 7335, NMAX1,DQSCA,DQEXT
583  IF(dsca.LE.ddelt.AND.dext.LE.ddelt) GO TO 55
584  IF (nma.EQ.npn1) print 7333, npn1
585  IF (nma.EQ.npn1) stop
586  50 CONTINUE
587  55 nnnggg=ngauss+1
588  IF (ngauss.EQ.npng1) print 7336
589  mmax=nmax1
590  DO 150 ngaus=nnnggg,npng1
591  ngauss=ngaus
592  nggg=2*ngauss
593  7336 FORMAT('WARNING: NGAUSS=NPNG1')
594  7337 FORMAT(' NG=',i3,' DC2=',d8.2,' DC1=',d8.2)
595  CALL const(ngauss,nmax,mmax,p,x,w,an,ann,s,ss,np,eps)
596  CALL vary(lam,mrr,mri,a,eps,np,ngauss,x,p,ppi,pir,pii,r,
597  & dr,ddr,drr,dri,nmax)
598  CALL tmatr0 (ngauss,x,w,an,ann,s,ss,ppi,pir,pii,r,dr,
599  & ddr,drr,dri,nmax,ncheck)
600  qext=0d0
601  qsca=0d0
602  DO 104 n=1,nmax
603  n1=n+nmax
604  tr1nn=tr1(n,n)
605  ti1nn=ti1(n,n)
606  tr1nn1=tr1(n1,n1)
607  ti1nn1=ti1(n1,n1)
608  dn1=dfloat(2*n+1)
609  qsca=qsca+dn1*(tr1nn*tr1nn+ti1nn*ti1nn
610  & +tr1nn1*tr1nn1+ti1nn1*ti1nn1)
611  qext=qext+(tr1nn+tr1nn1)*dn1
612  104 CONTINUE
613  dsca=abs((qsca1-qsca)/qsca)
614  dext=abs((qext1-qext)/qext)
615 c PRINT 7337, NGGG,DSCA,DEXT
616  qext1=qext
617  qsca1=qsca
618  IF(dsca.LE.ddelt.AND.dext.LE.ddelt) GO TO 155
619  IF (ngaus.EQ.npng1) print 7336
620  150 CONTINUE
621  155 CONTINUE
622  qsca=0d0
623  qext=0d0
624  nnm=nmax*2
625  DO 204 n=1,nnm
626  qext=qext+tr1(n,n)
627  204 CONTINUE
628  IF (nmax1.GT.npn4) print 7550, nmax1
629  7550 FORMAT ('NMAX1 = ',i3, ', i.e. greater than NPN4.',
630  & ' Execution terminated')
631  IF (nmax1.GT.npn4) stop
632  DO 213 n2=1,nmax1
633  nn2=n2+nmax
634  DO 213 n1=1,nmax1
635  nn1=n1+nmax
636  zz1=tr1(n1,n2)
637  rt11(1,n1,n2)=zz1
638  zz2=ti1(n1,n2)
639  it11(1,n1,n2)=zz2
640  zz3=tr1(n1,nn2)
641  rt12(1,n1,n2)=zz3
642  zz4=ti1(n1,nn2)
643  it12(1,n1,n2)=zz4
644  zz5=tr1(nn1,n2)
645  rt21(1,n1,n2)=zz5
646  zz6=ti1(nn1,n2)
647  it21(1,n1,n2)=zz6
648  zz7=tr1(nn1,nn2)
649  rt22(1,n1,n2)=zz7
650  zz8=ti1(nn1,nn2)
651  it22(1,n1,n2)=zz8
652  qsca=qsca+zz1*zz1+zz2*zz2+zz3*zz3+zz4*zz4
653  & +zz5*zz5+zz6*zz6+zz7*zz7+zz8*zz8
654  213 CONTINUE
655 C PRINT 7800,0,ABS(QEXT),QSCA,NMAX
656  DO 220 m=1,nmax1
657  CALL tmatr(m,ngauss,x,w,an,ann,s,ss,ppi,pir,pii,r,dr,
658  & ddr,drr,dri,nmax,ncheck)
659  nm=nmax-m+1
660  nm1=nmax1-m+1
661  m1=m+1
662  qsc=0d0
663  DO 214 n2=1,nm1
664  nn2=n2+m-1
665  n22=n2+nm
666  DO 214 n1=1,nm1
667  nn1=n1+m-1
668  n11=n1+nm
669  zz1=tr1(n1,n2)
670  rt11(m1,nn1,nn2)=zz1
671  zz2=ti1(n1,n2)
672  it11(m1,nn1,nn2)=zz2
673  zz3=tr1(n1,n22)
674  rt12(m1,nn1,nn2)=zz3
675  zz4=ti1(n1,n22)
676  it12(m1,nn1,nn2)=zz4
677  zz5=tr1(n11,n2)
678  rt21(m1,nn1,nn2)=zz5
679  zz6=ti1(n11,n2)
680  it21(m1,nn1,nn2)=zz6
681  zz7=tr1(n11,n22)
682  rt22(m1,nn1,nn2)=zz7
683  zz8=ti1(n11,n22)
684  it22(m1,nn1,nn2)=zz8
685  qsc=qsc+(zz1*zz1+zz2*zz2+zz3*zz3+zz4*zz4
686  & +zz5*zz5+zz6*zz6+zz7*zz7+zz8*zz8)*2d0
687  214 CONTINUE
688  nnm=2*nm
689  qxt=0d0
690  DO 215 n=1,nnm
691  qxt=qxt+tr1(n,n)*2d0
692  215 CONTINUE
693  qsca=qsca+qsc
694  qext=qext+qxt
695 C PRINT 7800,M,ABS(QXT),QSC,NMAX
696  7800 FORMAT(' m=',i3,' qxt=',d12.6,' qsc=',d12.6,
697  & ' nmax=',i3)
698  220 CONTINUE
699  coeff1=lam*lam*0.5d0/p
700  csca=qsca*coeff1
701  cext=-qext*coeff1
702 c PRINT 7880, NMAX,NMAX1
703  7880 FORMAT ('nmax=',i3,' nmax1=',i3)
704  CALL gsp (nmax1,csca,lam,al1,al2,al3,al4,be1,be2,lmax)
705  l1m=lmax+1
706  l1max=max(l1max,l1m)
707  wgii=wg1(i)
708  wgi=wgii*csca
709  DO 250 l1=1,l1m
710  alph1(l1)=alph1(l1)+al1(l1)*wgi
711  alph2(l1)=alph2(l1)+al2(l1)*wgi
712  alph3(l1)=alph3(l1)+al3(l1)*wgi
713  alph4(l1)=alph4(l1)+al4(l1)*wgi
714  bet1(l1)=bet1(l1)+be1(l1)*wgi
715  bet2(l1)=bet2(l1)+be2(l1)*wgi
716  250 CONTINUE
717  cscat=cscat+wgi
718  cextin=cextin+cext*wgii
719 C PRINT 6070, I,NMAX,NMAX1,NGAUSS
720  6070 FORMAT(4i6)
721  500 CONTINUE
722  DO 510 l1=1,l1max
723  alph1(l1)=alph1(l1)/cscat
724  alph2(l1)=alph2(l1)/cscat
725  alph3(l1)=alph3(l1)/cscat
726  alph4(l1)=alph4(l1)/cscat
727  bet1(l1)=bet1(l1)/cscat
728  bet2(l1)=bet2(l1)/cscat
729  alpha(l1,1)=alph1(l1)
730  alpha(l1,2)=alph2(l1)
731  alpha(l1,3)=alph3(l1)
732  alpha(l1,4)=alph4(l1)
733  beta(l1,1)=bet1(l1)
734  beta(l1,2)=bet2(l1)
735  510 CONTINUE
736  walb=cscat/cextin
737  CALL hovenr(l1max,alph1,alph2,alph3,alph4,bet1,bet2)
738  asymm=alph1(2)/3d0
739  print 9100,cextin,cscat,walb,asymm
740  9100 FORMAT('CEXT=',d12.6,2x,'CSCA=',d12.6,2x,
741  & 2x,'W=',d12.6,2x,'<COS>=',d12.6)
742  IF (walb.GT.1d0) print 9111
743  9111 FORMAT ('WARNING: W IS GREATER THAN 1')
744  WRITE (10,580) walb,l1max
745  DO l=1,l1max
746  WRITE (10,575) alph1(l),alph2(l),alph3(l),alph4(l),
747  & bet1(l),bet2(l)
748  ENDDO
749  575 FORMAT(6d14.7)
750  580 FORMAT(d14.8,i8)
751  lmax=l1max-1
752  CALL matr (alph1,alph2,alph3,alph4,bet1,bet2,lmax,npna,f)
753  600 CONTINUE
754  itime=mclock()
755  time=dfloat(itime)/6000d0
756  print 1001,time
757  1001 FORMAT (' time =',f8.2,' min')
758  RETURN
759  END
760 
761 C**********************************************************************
762 
763 C**********************************************************************
764 
765 C**********************************************************************
766 
767 C**********************************************************************
768 
769 C**********************************************************************
770 
771 C**********************************************************************
772 
773 C**********************************************************************
774 
775 C**********************************************************************
776 
777 C**********************************************************************
778 C *
779 C CALCULATION OF SPHERICAL BESSEL FUNCTIONS OF THE FIRST KIND *
780 C J=JR+I*JI OF COMPLEX ARGUMENT X=XR+I*XI OF ORDERS FROM 1 TO NMAX *
781 C BY USING BACKWARD RECURSION. PARAMETR NNMAX DETERMINES NUMERICAL *
782 C ACCURACY. U=UR+I*UI - FUNCTION (1/X)(D/DX)(X*J(X)) *
783 C *
784 C**********************************************************************
785 
786 C**********************************************************************
787 
788 C**********************************************************************
789 
790 C**********************************************************************
791 
792 C********************************************************************
793 C *
794 C CALCULATION OF THE EXPANSION COEFFICIENTS FOR (I,Q,U,V) - *
795 C REPRESENTATION. *
796 C *
797 C INPUT PARAMETERS: *
798 C *
799 C LAM - WAVELENGTH OF LIGHT *
800 C CSCA - SCATTERING CROSS SECTION *
801 C TR AND TI - ELEMENTS OF THE T-MATRIX. TRANSFERRED THROUGH *
802 C COMMON /CTM/ *
803 C NMAX - DIMENSION OF T(M)-MATRICES *
804 C *
805 C OUTPUT INFORTMATION: *
806 C *
807 C ALF1,...,ALF4,BET1,BET2 - EXPANSION COEFFICIENTS *
808 C LMAX - NUMBER OF COEFFICIENTS MINUS 1 *
809 C *
810 C********************************************************************
811 
812  SUBROUTINE gsp(NMAX,CSCA,LAM,ALF1,ALF2,ALF3,ALF4,BET1,BET2,LMAX)
813  include 'tmd.par.f'
814  IMPLICIT REAL*8 (a-b,d-h,o-z),COMPLEX*16 (C)
815  real*8 lam,ssign(900)
816  real*8 csca,ssi(npl),ssj(npn1),
817  & alf1(npl),alf2(npl),alf3(npl),
818  & alf4(npl),bet1(npl),bet2(npl),
819  & tr1(npl1,npn4),tr2(npl1,npn4),
820  & ti1(npl1,npn4),ti2(npl1,npn4),
821  & g1(npl1,npn6),g2(npl1,npn6),
822  & ar1(npn4),ar2(npn4),ai1(npn4),ai2(npn4),
823  & fr(npn4,npn4),fi(npn4,npn4),ff(npn4,npn4)
824  real*4 b1r(npl1,npl1,npn4),b1i(npl1,npl1,npn4),
825  & b2r(npl1,npl1,npn4),b2i(npl1,npl1,npn4),
826  & d1(npl1,npn4,npn4),d2(npl1,npn4,npn4),
827  & d3(npl1,npn4,npn4),d4(npl1,npn4,npn4),
828  & d5r(npl1,npn4,npn4),d5i(npl1,npn4,npn4),
829  & plus1(npn6*npn4*npn4*8)
830  real*4
831  & tr11(npn6,npn4,npn4),tr12(npn6,npn4,npn4),
832  & tr21(npn6,npn4,npn4),tr22(npn6,npn4,npn4),
833  & ti11(npn6,npn4,npn4),ti12(npn6,npn4,npn4),
834  & ti21(npn6,npn4,npn4),ti22(npn6,npn4,npn4)
835  COMPLEX*16 CIM(NPN1)
836 
837  COMMON /tmat/ tr11,tr12,tr21,tr22,ti11,ti12,ti21,ti22
838  COMMON /cbess/ b1r,b1i,b2r,b2i
839  COMMON /ss/ ssign
840  equivalence( plus1(1),tr11(1,1,1) )
841  equivalence(d1(1,1,1),plus1(1)),
842  & (d2(1,1,1),plus1(npl1*npn4*npn4+1)),
843  & (d3(1,1,1),plus1(npl1*npn4*npn4*2+1)),
844  & (d4(1,1,1),plus1(npl1*npn4*npn4*3+1)),
845  & (d5r(1,1,1),plus1(npl1*npn4*npn4*4+1))
846 
847  CALL fact
848  CALL signum
849  lmax=2*nmax
850  l1max=lmax+1
851  ci=(0d0,1d0)
852  cim(1)=ci
853  DO 2 i=2,nmax
854  cim(i)=cim(i-1)*ci
855  2 CONTINUE
856  ssi(1)=1d0
857  DO 3 i=1,lmax
858  i1=i+1
859  si=dfloat(2*i+1)
860  ssi(i1)=si
861  IF(i.LE.nmax) ssj(i)=dsqrt(si)
862  3 CONTINUE
863  ci=-ci
864  DO 5 i=1,nmax
865  si=ssj(i)
866  cci=cim(i)
867  DO 4 j=1,nmax
868  sj=1d0/ssj(j)
869  ccj=cim(j)*sj/cci
870  fr(j,i)=ccj
871  fi(j,i)=ccj*ci
872  ff(j,i)=si*sj
873  4 CONTINUE
874  5 CONTINUE
875  nmax1=nmax+1
876 
877 C ***** CALCULATION OF THE ARRAYS B1 AND B2 *****
878 
879  k1=1
880  k2=0
881  k3=0
882  k4=1
883  k5=1
884  k6=2
885 
886 C PRINT 3300, B1,B2
887  3300 FORMAT (' B1 AND B2')
888  DO 100 n=1,nmax
889 
890 C ***** CALCULATION OF THE ARRAYS T1 AND T2 *****
891 
892 
893  DO 10 nn=1,nmax
894  m1max=min0(n,nn)+1
895  DO 6 m1=1,m1max
896  m=m1-1
897  l1=npn6+m
898  tt1=tr11(m1,n,nn)
899  tt2=tr12(m1,n,nn)
900  tt3=tr21(m1,n,nn)
901  tt4=tr22(m1,n,nn)
902  tt5=ti11(m1,n,nn)
903  tt6=ti12(m1,n,nn)
904  tt7=ti21(m1,n,nn)
905  tt8=ti22(m1,n,nn)
906  t1=tt1+tt2
907  t2=tt3+tt4
908  t3=tt5+tt6
909  t4=tt7+tt8
910  tr1(l1,nn)=t1+t2
911  tr2(l1,nn)=t1-t2
912  ti1(l1,nn)=t3+t4
913  ti2(l1,nn)=t3-t4
914  IF(m.EQ.0) GO TO 6
915  l1=npn6-m
916  t1=tt1-tt2
917  t2=tt3-tt4
918  t3=tt5-tt6
919  t4=tt7-tt8
920  tr1(l1,nn)=t1-t2
921  tr2(l1,nn)=t1+t2
922  ti1(l1,nn)=t3-t4
923  ti2(l1,nn)=t3+t4
924  6 CONTINUE
925  10 CONTINUE
926 
927 C ***** END OF THE CALCULATION OF THE ARRAYS T1 AND T2 *****
928 
929  nn1max=nmax1+n
930  DO 40 nn1=1,nn1max
931  n1=nn1-1
932 
933 C ***** CALCULATION OF THE ARRAYS A1 AND A2 *****
934 
935  CALL ccg(n,n1,nmax,k1,k2,g1)
936  nnmax=min0(nmax,n1+n)
937  nnmin=max0(1,iabs(n-n1))
938  kn=n+nn1
939  DO 15 nn=nnmin,nnmax
940  nnn=nn+1
941  sig=ssign(kn+nn)
942  m1max=min0(n,nn)+npn6
943  aar1=0d0
944  aar2=0d0
945  aai1=0d0
946  aai2=0d0
947  DO 13 m1=npn6,m1max
948  m=m1-npn6
949  sss=g1(m1,nnn)
950  rr1=tr1(m1,nn)
951  ri1=ti1(m1,nn)
952  rr2=tr2(m1,nn)
953  ri2=ti2(m1,nn)
954  IF(m.EQ.0) GO TO 12
955  m2=npn6-m
956  rr1=rr1+tr1(m2,nn)*sig
957  ri1=ri1+ti1(m2,nn)*sig
958  rr2=rr2+tr2(m2,nn)*sig
959  ri2=ri2+ti2(m2,nn)*sig
960  12 aar1=aar1+sss*rr1
961  aai1=aai1+sss*ri1
962  aar2=aar2+sss*rr2
963  aai2=aai2+sss*ri2
964  13 CONTINUE
965  xr=fr(nn,n)
966  xi=fi(nn,n)
967  ar1(nn)=aar1*xr-aai1*xi
968  ai1(nn)=aar1*xi+aai1*xr
969  ar2(nn)=aar2*xr-aai2*xi
970  ai2(nn)=aar2*xi+aai2*xr
971  15 CONTINUE
972 
973 C ***** END OF THE CALCULATION OF THE ARRAYS A1 AND A2 ****
974 
975  CALL ccg(n,n1,nmax,k3,k4,g2)
976  m1=max0(-n1+1,-n)
977  m2=min0(n1+1,n)
978  m1max=m2+npn6
979  m1min=m1+npn6
980  DO 30 m1=m1min,m1max
981  bbr1=0d0
982  bbi1=0d0
983  bbr2=0d0
984  bbi2=0d0
985  DO 25 nn=nnmin,nnmax
986  nnn=nn+1
987  sss=g2(m1,nnn)
988  bbr1=bbr1+sss*ar1(nn)
989  bbi1=bbi1+sss*ai1(nn)
990  bbr2=bbr2+sss*ar2(nn)
991  bbi2=bbi2+sss*ai2(nn)
992  25 CONTINUE
993  b1r(nn1,m1,n)=bbr1
994  b1i(nn1,m1,n)=bbi1
995  b2r(nn1,m1,n)=bbr2
996  b2i(nn1,m1,n)=bbi2
997  30 CONTINUE
998  40 CONTINUE
999  100 CONTINUE
1000 
1001 C ***** END OF THE CALCULATION OF THE ARRAYS B1 AND B2 ****
1002 
1003 C ***** CALCULATION OF THE ARRAYS D1,D2,D3,D4, AND D5 *****
1004 
1005 C PRINT 3301
1006  3301 FORMAT(' D1, D2, ...')
1007  DO 200 n=1,nmax
1008  DO 190 nn=1,nmax
1009  m1=min0(n,nn)
1010  m1max=npn6+m1
1011  m1min=npn6-m1
1012  nn1max=nmax1+min0(n,nn)
1013  DO 180 m1=m1min,m1max
1014  m=m1-npn6
1015  nn1min=iabs(m-1)+1
1016  dd1=0d0
1017  dd2=0d0
1018  DO 150 nn1=nn1min,nn1max
1019  xx=ssi(nn1)
1020  x1=b1r(nn1,m1,n)
1021  x2=b1i(nn1,m1,n)
1022  x3=b1r(nn1,m1,nn)
1023  x4=b1i(nn1,m1,nn)
1024  x5=b2r(nn1,m1,n)
1025  x6=b2i(nn1,m1,n)
1026  x7=b2r(nn1,m1,nn)
1027  x8=b2i(nn1,m1,nn)
1028  dd1=dd1+xx*(x1*x3+x2*x4)
1029  dd2=dd2+xx*(x5*x7+x6*x8)
1030  150 CONTINUE
1031  d1(m1,nn,n)=dd1
1032  d2(m1,nn,n)=dd2
1033  180 CONTINUE
1034  mmax=min0(n,nn+2)
1035  mmin=max0(-n,-nn+2)
1036  m1max=npn6+mmax
1037  m1min=npn6+mmin
1038  DO 186 m1=m1min,m1max
1039  m=m1-npn6
1040  nn1min=iabs(m-1)+1
1041  dd3=0d0
1042  dd4=0d0
1043  dd5r=0d0
1044  dd5i=0d0
1045  m2=-m+2+npn6
1046  DO 183 nn1=nn1min,nn1max
1047  xx=ssi(nn1)
1048  x1=b1r(nn1,m1,n)
1049  x2=b1i(nn1,m1,n)
1050  x3=b2r(nn1,m1,n)
1051  x4=b2i(nn1,m1,n)
1052  x5=b1r(nn1,m2,nn)
1053  x6=b1i(nn1,m2,nn)
1054  x7=b2r(nn1,m2,nn)
1055  x8=b2i(nn1,m2,nn)
1056  dd3=dd3+xx*(x1*x5+x2*x6)
1057  dd4=dd4+xx*(x3*x7+x4*x8)
1058  dd5r=dd5r+xx*(x3*x5+x4*x6)
1059  dd5i=dd5i+xx*(x4*x5-x3*x6)
1060  183 CONTINUE
1061  d3(m1,nn,n)=dd3
1062  d4(m1,nn,n)=dd4
1063  d5r(m1,nn,n)=dd5r
1064  d5i(m1,nn,n)=dd5i
1065  186 CONTINUE
1066  190 CONTINUE
1067  200 CONTINUE
1068 
1069 C ***** END OF THE CALCULATION OF THE D-ARRAYS *****
1070 
1071 C ***** CALCULATION OF THE EXPANSION COEFFICIENTS *****
1072 
1073 C PRINT 3303
1074  3303 FORMAT (' G1, G2, ...')
1075 
1076  dk=lam*lam/(4d0*csca*dacos(-1d0))
1077  DO 300 l1=1,l1max
1078  g1l=0d0
1079  g2l=0d0
1080  g3l=0d0
1081  g4l=0d0
1082  g5lr=0d0
1083  g5li=0d0
1084  l=l1-1
1085  sl=ssi(l1)*dk
1086  DO 290 n=1,nmax
1087  nnmin=max0(1,iabs(n-l))
1088  nnmax=min0(nmax,n+l)
1089  IF(nnmax.LT.nnmin) GO TO 290
1090  CALL ccg(n,l,nmax,k1,k2,g1)
1091  IF(l.GE.2) CALL ccg(n,l,nmax,k5,k6,g2)
1092  nl=n+l
1093  DO 280 nn=nnmin,nnmax
1094  nnn=nn+1
1095  mmax=min0(n,nn)
1096  m1min=npn6-mmax
1097  m1max=npn6+mmax
1098  si=ssign(nl+nnn)
1099  dm1=0d0
1100  dm2=0d0
1101  DO 270 m1=m1min,m1max
1102  m=m1-npn6
1103  IF(m.GE.0) sss1=g1(m1,nnn)
1104  IF(m.LT.0) sss1=g1(npn6-m,nnn)*si
1105  dm1=dm1+sss1*d1(m1,nn,n)
1106  dm2=dm2+sss1*d2(m1,nn,n)
1107  270 CONTINUE
1108  ffn=ff(nn,n)
1109  sss=g1(npn6+1,nnn)*ffn
1110  g1l=g1l+sss*dm1
1111  g2l=g2l+sss*dm2*si
1112  IF(l.LT.2) GO TO 280
1113  dm3=0d0
1114  dm4=0d0
1115  dm5r=0d0
1116  dm5i=0d0
1117  mmax=min0(n,nn+2)
1118  mmin=max0(-n,-nn+2)
1119  m1max=npn6+mmax
1120  m1min=npn6+mmin
1121  DO 275 m1=m1min,m1max
1122  m=m1-npn6
1123  sss1=g2(npn6-m,nnn)
1124  dm3=dm3+sss1*d3(m1,nn,n)
1125  dm4=dm4+sss1*d4(m1,nn,n)
1126  dm5r=dm5r+sss1*d5r(m1,nn,n)
1127  dm5i=dm5i+sss1*d5i(m1,nn,n)
1128  275 CONTINUE
1129  g5lr=g5lr-sss*dm5r
1130  g5li=g5li-sss*dm5i
1131  sss=g2(npn4,nnn)*ffn
1132  g3l=g3l+sss*dm3
1133  g4l=g4l+sss*dm4*si
1134  280 CONTINUE
1135  290 CONTINUE
1136  g1l=g1l*sl
1137  g2l=g2l*sl
1138  g3l=g3l*sl
1139  g4l=g4l*sl
1140  g5lr=g5lr*sl
1141  g5li=g5li*sl
1142  alf1(l1)=g1l+g2l
1143  alf2(l1)=g3l+g4l
1144  alf3(l1)=g3l-g4l
1145  alf4(l1)=g1l-g2l
1146  bet1(l1)=g5lr*2d0
1147  bet2(l1)=g5li*2d0
1148  lmax=l
1149  IF(abs(g1l).LT.1d-6) GO TO 500
1150  300 CONTINUE
1151  500 CONTINUE
1152  RETURN
1153  END
1154 
1155 C****************************************************************
1156 
1157 C CALCULATION OF THE QUANTITIES F(N+1)=0.5*LN(N!)
1158 C 0.LE.N.LE.899
1159 
1160  SUBROUTINE fact
1161  real*8 f(900)
1162  COMMON /fac/ f
1163  f(1)=0d0
1164  f(2)=0d0
1165  DO 2 i=3,900
1166  i1=i-1
1167  f(i)=f(i1)+0.5d0*dlog(dfloat(i1))
1168  2 CONTINUE
1169  RETURN
1170  END
1171 
1172 C************************************************************
1173 
1174 C CALCULATION OF THE ARRAY SSIGN(N+1)=SIGN(N)
1175 C 0.LE.N.LE.899
1176 
1177  SUBROUTINE signum
1178  real*8 ssign(900)
1179  COMMON /ss/ ssign
1180  ssign(1)=1d0
1181  DO 2 n=2,899
1182  ssign(n)=-ssign(n-1)
1183  2 CONTINUE
1184  RETURN
1185  END
1186 
1187 C******************************************************************
1188 C
1189 C CALCULATION OF CLEBSCH-GORDAN COEFFICIENTS
1190 C (N,M:N1,M1/NN,MM)
1191 C FOR GIVEN N AND N1. M1=MM-M, INDEX MM IS FOUND FROM M AS
1192 C MM=M*K1+K2
1193 C
1194 C INPUT PARAMETERS : N,N1,NMAX,K1,K2
1195 C N.LE.NMAX
1196 C N.GE.1
1197 C N1.GE.0
1198 C N1.LE.N+NMAX
1199 C OUTPUT PARAMETERS : GG(M+NPN6,NN+1) - ARRAY OF THE CORRESPONDING
1200 C COEFFICIENTS
1201 C /M/.LE.N
1202 C /M1/=/M*(K1-1)+K2/.LE.N1
1203 C NN.LE.MIN(N+N1,NMAX)
1204 C NN.GE.MAX(/MM/,/N-N1/)
1205 C IF K1=1 AND K2=0, THEN 0.LE.M.LE.N
1206 
1207 
1208  SUBROUTINE ccg(N,N1,NMAX,K1,K2,GG)
1209  include 'tmd.par.f'
1210  IMPLICIT REAL*8 (a-h,o-z)
1211  real*8 gg(npl1,npn6),cd(0:npn5),cu(0:npn5)
1212  IF(nmax.LE.npn4.
1213  & and.0.LE.n1.
1214  & and.n1.LE.nmax+n.
1215  & and.n.GE.1.
1216  & and.n.LE.nmax) GO TO 1
1217  print 5001
1218  stop
1219  5001 FORMAT(' ERROR IN SUBROUTINE CCG')
1220  1 nnf=min0(n+n1,nmax)
1221  min=npn6-n
1222  mf=npn6+n
1223  IF(k1.EQ.1.AND.k2.EQ.0) min=npn6
1224  DO 100 mind=min,mf
1225  m=mind-npn6
1226  mm=m*k1+k2
1227  m1=mm-m
1228  IF(iabs(m1).GT.n1) GO TO 90
1229  nnl=max0(iabs(mm),iabs(n-n1))
1230  IF(nnl.GT.nnf) GO TO 90
1231  nnu=n+n1
1232  nnm=(nnu+nnl)*0.5d0
1233  IF (nnu.EQ.nnl) nnm=nnl
1234  CALL ccgin(n,n1,m,mm,c)
1235  cu(nnl)=c
1236  IF (nnl.EQ.nnf) GO TO 50
1237  c2=0d0
1238  c1=c
1239  DO 7 nn=nnl+1,min0(nnm,nnf)
1240  a=dfloat((nn+mm)*(nn-mm)*(n1-n+nn))
1241  a=a*dfloat((n-n1+nn)*(n+n1-nn+1)*(n+n1+nn+1))
1242  a=dfloat(4*nn*nn)/a
1243  a=a*dfloat((2*nn+1)*(2*nn-1))
1244  a=dsqrt(a)
1245  b=0.5d0*dfloat(m-m1)
1246  d=0d0
1247  IF(nn.EQ.1) GO TO 5
1248  b=dfloat(2*nn*(nn-1))
1249  b=dfloat((2*m-mm)*nn*(nn-1)-mm*n*(n+1)+
1250  & mm*n1*(n1+1))/b
1251  d=dfloat(4*(nn-1)*(nn-1))
1252  d=d*dfloat((2*nn-3)*(2*nn-1))
1253  d=dfloat((nn-mm-1)*(nn+mm-1)*(n1-n+nn-1))/d
1254  d=d*dfloat((n-n1+nn-1)*(n+n1-nn+2)*(n+n1+nn))
1255  d=dsqrt(d)
1256  5 c=a*(b*c1-d*c2)
1257  c2=c1
1258  c1=c
1259  cu(nn)=c
1260  7 CONTINUE
1261  IF (nnf.LE.nnm) GO TO 50
1262  CALL direct(n,m,n1,m1,nnu,mm,c)
1263  cd(nnu)=c
1264  IF (nnu.EQ.nnm+1) GO TO 50
1265  c2=0d0
1266  c1=c
1267  DO 12 nn=nnu-1,nnm+1,-1
1268  a=dfloat((nn-mm+1)*(nn+mm+1)*(n1-n+nn+1))
1269  a=a*dfloat((n-n1+nn+1)*(n+n1-nn)*(n+n1+nn+2))
1270  a=dfloat(4*(nn+1)*(nn+1))/a
1271  a=a*dfloat((2*nn+1)*(2*nn+3))
1272  a=dsqrt(a)
1273  b=dfloat(2*(nn+2)*(nn+1))
1274  b=dfloat((2*m-mm)*(nn+2)*(nn+1)-mm*n*(n+1)
1275  & +mm*n1*(n1+1))/b
1276  d=dfloat(4*(nn+2)*(nn+2))
1277  d=d*dfloat((2*nn+5)*(2*nn+3))
1278  d=dfloat((nn+mm+2)*(nn-mm+2)*(n1-n+nn+2))/d
1279  d=d*dfloat((n-n1+nn+2)*(n+n1-nn-1)*(n+n1+nn+3))
1280  d=dsqrt(d)
1281  c=a*(b*c1-d*c2)
1282  c2=c1
1283  c1=c
1284  cd(nn)=c
1285  12 CONTINUE
1286  50 DO 9 nn=nnl,nnf
1287  IF (nn.LE.nnm) gg(mind,nn+1)=cu(nn)
1288  IF (nn.GT.nnm) gg(mind,nn+1)=cd(nn)
1289 c WRITE (6,*) N,M,N1,M1,NN,MM,GG(MIND,NN+1)
1290  9 CONTINUE
1291  90 CONTINUE
1292  100 CONTINUE
1293  RETURN
1294  END
1295 
1296 C*********************************************************************
1297 
1298  SUBROUTINE direct (N,M,N1,M1,NN,MM,C)
1299  IMPLICIT REAL*8 (a-h,o-z)
1300  real*8 f(900)
1301  COMMON /fac/ f
1302  c=f(2*n+1)+f(2*n1+1)+f(n+n1+m+m1+1)+f(n+n1-m-m1+1)
1303  c=c-f(2*(n+n1)+1)-f(n+m+1)-f(n-m+1)-f(n1+m1+1)-f(n1-m1+1)
1304  c=dexp(c)
1305  RETURN
1306  END
1307 
1308 C*********************************************************************
1309 C
1310 C CALCULATION OF THE CLEBCSH-GORDAN COEFFICIENTS
1311 C G=(N,M:N1,MM-M/NN,MM)
1312 C FOR GIVEN N,N1,M,MM, WHERE NN=MAX(/MM/,/N-N1/)
1313 C /M/.LE.N
1314 C /MM-M/.LE.N1
1315 C /MM/.LE.N+N1
1316 
1317  SUBROUTINE ccgin(N,N1,M,MM,G)
1318  IMPLICIT REAL*8 (a-h,o-z)
1319  real*8 f(900),ssign(900)
1320  COMMON /ss/ ssign
1321  COMMON /fac/ f
1322  m1=mm-m
1323  IF(n.GE.iabs(m).
1324  & and.n1.GE.iabs(m1).
1325  & and.iabs(mm).LE.(n+n1)) GO TO 1
1326  print 5001
1327  stop
1328  5001 FORMAT(' ERROR IN SUBROUTINE CCGIN')
1329  1 IF (iabs(mm).GT.iabs(n-n1)) GO TO 100
1330  l1=n
1331  l2=n1
1332  l3=m
1333  IF(n1.LE.n) GO TO 50
1334  k=n
1335  n=n1
1336  n1=k
1337  k=m
1338  m=m1
1339  m1=k
1340  50 n2=n*2
1341  m2=m*2
1342  n12=n1*2
1343  m12=m1*2
1344  g=ssign(n1+m1+1)
1345  & *dexp(f(n+m+1)+f(n-m+1)+f(n12+1)+f(n2-n12+2)-f(n2+2)
1346  & -f(n1+m1+1)-f(n1-m1+1)-f(n-n1+mm+1)-f(n-n1-mm+1))
1347  n=l1
1348  n1=l2
1349  m=l3
1350  RETURN
1351  100 a=1d0
1352  l1=m
1353  l2=mm
1354  IF(mm.GE.0) GO TO 150
1355  mm=-mm
1356  m=-m
1357  m1=-m1
1358  a=ssign(mm+n+n1+1)
1359  150 g=a*ssign(n+m+1)
1360  & *dexp(f(2*mm+2)+f(n+n1-mm+1)+f(n+m+1)+f(n1+m1+1)
1361  & -f(n+n1+mm+2)-f(n-n1+mm+1)-f(-n+n1+mm+1)-f(n-m+1)
1362  & -f(n1-m1+1))
1363  m=l1
1364  mm=l2
1365  RETURN
1366  END
1367 
1368 C********************************************************************
1369 
1370 c********************************************************************
1371 
1372 C********************************************************************
1373 
1374 C********************************************************************
1375 
1376 C COMPUTATION OF R1 AND R2 FOR A POWER LAW SIZE DISTRIBUTION WITH
1377 C EFFECTIVE RADIUS A AND EFFECTIVE VARIANCE B
1378 
1379  SUBROUTINE power (A,B,R1,R2)
1380  IMPLICIT REAL*8 (a-h,o-z)
1381  EXTERNAL f
1382  COMMON aa,bb
1383  aa=a
1384  bb=b
1385  ax=1d-5
1386  bx=a-1d-5
1387  r1=zeroin(ax,bx,f,0d0)
1388  r2=(1d0+b)*2d0*a-r1
1389  RETURN
1390  END
1391 
1392 C***********************************************************************
1393 
1394  DOUBLE PRECISION FUNCTION zeroin (AX,BX,F,TOL)
1395  IMPLICIT REAL*8 (a-h,o-z)
1396  eps=1d0
1397  10 eps=0.5d0*eps
1398  tol1=1d0+eps
1399  IF (tol1.GT.1d0) GO TO 10
1400  15 a=ax
1401  b=bx
1402  fa=f(a)
1403  fb=f(b)
1404  20 c=a
1405  fc=fa
1406  d=b-a
1407  e=d
1408  30 IF (abs(fc).GE.abs(fb)) GO TO 40
1409  35 a=b
1410  b=c
1411  c=a
1412  fa=fb
1413  fb=fc
1414  fc=fa
1415  40 tol1=2d0*eps*abs(b)+0.5d0*tol
1416  xm=0.5d0*(c-b)
1417  IF (abs(xm).LE.tol1) GO TO 90
1418  44 IF (fb.EQ.0d0) GO TO 90
1419  45 IF (abs(e).LT.tol1) GO TO 70
1420  46 IF (abs(fa).LE.abs(fb)) GO TO 70
1421  47 IF (a.NE.c) GO TO 50
1422  48 s=fb/fa
1423  p=2d0*xm*s
1424  q=1d0-s
1425  GO TO 60
1426  50 q=fa/fc
1427  r=fb/fc
1428  s=fb/fa
1429  p=s*(2d0*xm*q*(q-r)-(b-a)*(r-1d0))
1430  q=(q-1d0)*(r-1d0)*(s-1d0)
1431  60 IF (p.GT.0d0) q=-q
1432  p=abs(p)
1433  IF ((2d0*p).GE.(3d0*xm*q-abs(tol1*q))) GO TO 70
1434  64 IF (p.GE.abs(0.5d0*e*q)) GO TO 70
1435  65 e=d
1436  d=p/q
1437  GO TO 80
1438  70 d=xm
1439  e=d
1440  80 a=b
1441  fa=fb
1442  IF (abs(d).GT.tol1) b=b+d
1443  IF (abs(d).LE.tol1) b=b+dsign(tol1,xm)
1444  fb=f(b)
1445  IF ((fb*(fc/abs(fc))).GT.0d0) GO TO 20
1446  85 GO TO 30
1447  90 zeroin=b
1448  RETURN
1449  END
1450 
1451 C***********************************************************************
1452 
1453  DOUBLE PRECISION FUNCTION f(R1)
1454  IMPLICIT REAL*8 (a-h,o-z)
1455  COMMON a,b
1456  r2=(1d0+b)*2d0*a-r1
1457  f=(r2-r1)/dlog(r2/r1)-a
1458  RETURN
1459  END
1460 
1461 C**********************************************************************
1462 C CALCULATION OF POINTS AND WEIGHTS OF GAUSSIAN QUADRATURE *
1463 C FORMULA. IF IND1 = 0 - ON INTERVAL (-1,1), IF IND1 = 1 - ON *
1464 C INTERVAL (0,1). IF IND2 = 1 RESULTS ARE PRINTED. *
1465 C N - NUMBER OF POINTS *
1466 C Z - DIVISION POINTS *
1467 C W - WEIGHTS *
1468 C**********************************************************************
1469 
1470 C****************************************************************
1471 
1472  SUBROUTINE distrb (NNK,YY,WY,NDISTR,AA,BB,GAM,R1,R2,REFF,
1473  & VEFF,PI)
1474  IMPLICIT REAL*8 (a-h,o-z)
1475  real*8 yy(nnk),wy(nnk)
1476  IF (ndistr.EQ.2) GO TO 100
1477  IF (ndistr.EQ.3) GO TO 200
1478  IF (ndistr.EQ.4) GO TO 300
1479  IF (ndistr.EQ.5) GO TO 360
1480  print 1001,aa,bb,gam
1481  1001 FORMAT('MODIFIED GAMMA DISTRIBUTION, alpha=',f6.4,' r_c=',
1482  & f6.4,' gamma=',f6.4)
1483  a2=aa/gam
1484  db=1d0/bb
1485  DO 50 i=1,nnk
1486  x=yy(i)
1487  y=x**aa
1488  x=x*db
1489  y=y*dexp(-a2*(x**gam))
1490  wy(i)=wy(i)*y
1491  50 CONTINUE
1492  GO TO 400
1493  100 print 1002,aa,bb
1494  1002 FORMAT('LOG-NORMAL DISTRIBUTION, r_g=',f8.4,
1495  & ' [ln(sigma_g)]**2=', f6.4)
1496  da=1d0/aa
1497  DO 150 i=1,nnk
1498  x=yy(i)
1499  y=dlog(x*da)
1500  y=dexp(-y*y*0.5d0/bb)/x
1501  wy(i)=wy(i)*y
1502  150 CONTINUE
1503  GO TO 400
1504  200 print 1003
1505  1003 FORMAT('POWER LAW DISTRIBUTION OF HANSEN & TRAVIS 1974')
1506  DO 250 i=1,nnk
1507  x=yy(i)
1508  wy(i)=wy(i)/(x*x*x)
1509  250 CONTINUE
1510  GO TO 400
1511  300 print 1004,aa,bb
1512  1004 FORMAT ('GAMMA DISTRIBUTION, a=',f6.3,' b=',f6.4)
1513  b2=(1d0-3d0*bb)/bb
1514  dab=1d0/(aa*bb)
1515  DO 350 i=1,nnk
1516  x=yy(i)
1517  x=(x**b2)*dexp(-x*dab)
1518  wy(i)=wy(i)*x
1519  350 CONTINUE
1520  GO TO 400
1521  360 print 1005,bb
1522  1005 FORMAT ('MODIFIED POWER LAW DISTRIBUTION, alpha=',d10.4)
1523  DO 370 i=1,nnk
1524  x=yy(i)
1525  IF (x.LE.r1) wy(i)=wy(i)
1526  IF (x.GT.r1) wy(i)=wy(i)*(x/r1)**bb
1527  370 CONTINUE
1528  400 CONTINUE
1529  sum=0d0
1530  DO 450 i=1,nnk
1531  sum=sum+wy(i)
1532  450 CONTINUE
1533  sum=1d0/sum
1534  DO 500 i=1,nnk
1535  wy(i)=wy(i)*sum
1536  500 CONTINUE
1537  g=0d0
1538  DO 550 i=1,nnk
1539  x=yy(i)
1540  g=g+x*x*wy(i)
1541  550 CONTINUE
1542  reff=0d0
1543  DO 600 i=1,nnk
1544  x=yy(i)
1545  reff=reff+x*x*x*wy(i)
1546  600 CONTINUE
1547  reff=reff/g
1548  veff=0d0
1549  DO 650 i=1,nnk
1550  x=yy(i)
1551  xi=x-reff
1552  veff=veff+xi*xi*x*x*wy(i)
1553  650 CONTINUE
1554  veff=veff/(g*reff*reff)
1555  RETURN
1556  END
1557 
1558 C*************************************************************
1559 
1560  SUBROUTINE hovenr(L1,A1,A2,A3,A4,B1,B2)
1561  IMPLICIT REAL*8 (a-h,o-z)
1562  real*8 a1(l1),a2(l1),a3(l1),a4(l1),b1(l1),b2(l1)
1563  DO 100 l=1,l1
1564  kontr=1
1565  ll=l-1
1566  dl=dfloat(ll)*2d0+1d0
1567  ddl=dl*0.48d0
1568  aa1=a1(l)
1569  aa2=a2(l)
1570  aa3=a3(l)
1571  aa4=a4(l)
1572  bb1=b1(l)
1573  bb2=b2(l)
1574  IF(ll.GE.1.AND.abs(aa1).GE.dl) kontr=2
1575  IF(abs(aa2).GE.dl) kontr=2
1576  IF(abs(aa3).GE.dl) kontr=2
1577  IF(abs(aa4).GE.dl) kontr=2
1578  IF(abs(bb1).GE.ddl) kontr=2
1579  IF(abs(bb2).GE.ddl) kontr=2
1580  IF(kontr.EQ.2) print 3000,ll
1581  c=-0.1d0
1582  DO 50 i=1,11
1583  c=c+0.1d0
1584  cc=c*c
1585  c1=cc*bb2*bb2
1586  c2=c*aa4
1587  c3=c*aa3
1588  IF((dl-c*aa1)*(dl-c*aa2)-cc*bb1*bb1.LE.-1d-4) kontr=2
1589  IF((dl-c2)*(dl-c3)+c1.LE.-1d-4) kontr=2
1590  IF((dl+c2)*(dl-c3)-c1.LE.-1d-4) kontr=2
1591  IF((dl-c2)*(dl+c3)-c1.LE.-1d-4) kontr=2
1592  IF(kontr.EQ.2) print 4000,ll,c
1593  50 CONTINUE
1594  100 CONTINUE
1595  IF(kontr.EQ.1) print 2000
1596  2000 FORMAT('TEST OF VAN DER MEE & HOVENIER IS SATISFIED')
1597  3000 FORMAT('TEST OF VAN DER MEE & HOVENIER IS NOT SATISFIED, L=',i3)
1598  4000 FORMAT('TEST OF VAN DER MEE & HOVENIER IS NOT SATISFIED, L=',i3,
1599  & ' A=',d9.2)
1600  RETURN
1601  END
1602 
1603 C****************************************************************
1604 
1605 C CALCULATION OF THE SCATTERING MATRIX FOR GIVEN EXPANSION
1606 C COEFFICIENTS
1607 
1608 C A1,...,B2 - EXPANSION COEFFICIENTS
1609 C LMAX - NUMBER OF COEFFICIENTS MINUS 1
1610 C N - NUMBER OF SCATTERING ANGLES
1611 C THE CORRESPONDING SCATTERING ANGLES ARE GIVEN BY
1612 C 180*(I-1)/(N-1) (DEGREES), WHERE I NUMBERS THE ANGLES
1613 
1614  SUBROUTINE matr(A1,A2,A3,A4,B1,B2,LMAX,NPNA,F)
1615  include 'tmd.par.f'
1616  IMPLICIT REAL*8 (a-h,o-z)
1617  real*8 a1(npl),a2(npl),a3(npl),a4(npl),b1(npl),b2(npl),f(npna,4,4)
1618  n=npna
1619  dn=1d0/dfloat(n-1)
1620  da=dacos(-1d0)*dn
1621  db=180d0*dn
1622  l1max=lmax+1
1623  print 1000
1624  1000 FORMAT(' ')
1625  print 1001
1626  1001 FORMAT(' ',2x,'S',6x,'ALPHA1',6x,'ALPHA2',6x,'ALPHA3',
1627  & 6x,'ALPHA4',7x,'BETA1',7x,'BETA2')
1628  DO 10 l1=1,l1max
1629  l=l1-1
1630  print 1002,l,a1(l1),a2(l1),a3(l1),a4(l1),b1(l1),b2(l1)
1631  10 CONTINUE
1632  1002 FORMAT(' ',i3,6f12.5)
1633  tb=-db
1634  taa=-da
1635  print 1000
1636  print 1003
1637  1003 FORMAT(' ',5x,'<',8x,'F11',8x,'F22',8x,'F33',
1638  & 8x,'F44',8x,'F12',8x,'F34')
1639  d6=dsqrt(6d0)*0.25d0
1640  DO 500 i1=1,n
1641  taa=taa+da
1642  tb=tb+db
1643  u=dcos(taa)
1644  f11=0d0
1645  f2=0d0
1646  f3=0d0
1647  f44=0d0
1648  f12=0d0
1649  f34=0d0
1650  p1=0d0
1651  p2=0d0
1652  p3=0d0
1653  p4=0d0
1654  pp1=1d0
1655  pp2=0.25d0*(1d0+u)*(1d0+u)
1656  pp3=0.25d0*(1d0-u)*(1d0-u)
1657  pp4=d6*(u*u-1d0)
1658  DO 400 l1=1,l1max
1659  l=l1-1
1660  dl=dfloat(l)
1661  dl1=dfloat(l1)
1662  f11=f11+a1(l1)*pp1
1663  f44=f44+a4(l1)*pp1
1664  IF(l.EQ.lmax) GO TO 350
1665  pl1=dfloat(2*l+1)
1666  p=(pl1*u*pp1-dl*p1)/dl1
1667  p1=pp1
1668  pp1=p
1669  350 IF(l.LT.2) GO TO 400
1670  f2=f2+(a2(l1)+a3(l1))*pp2
1671  f3=f3+(a2(l1)-a3(l1))*pp3
1672  f12=f12+b1(l1)*pp4
1673  f34=f34+b2(l1)*pp4
1674  IF(l.EQ.lmax) GO TO 400
1675  pl2=dfloat(l*l1)*u
1676  pl3=dfloat(l1*(l*l-4))
1677  pl4=1d0/dfloat(l*(l1*l1-4))
1678  p=(pl1*(pl2-4d0)*pp2-pl3*p2)*pl4
1679  p2=pp2
1680  pp2=p
1681  p=(pl1*(pl2+4d0)*pp3-pl3*p3)*pl4
1682  p3=pp3
1683  pp3=p
1684  p=(pl1*u*pp4-dsqrt(dfloat(l*l-4))*p4)/dsqrt(dfloat(l1*l1-4))
1685  p4=pp4
1686  pp4=p
1687  400 CONTINUE
1688  f22=(f2+f3)*0.5d0
1689  f33=(f2-f3)*0.5d0
1690 C F22=F22/F11
1691 C F33=F33/F11
1692 C F44=F44/F11
1693 C F12=-F12/F11
1694 C F34=F34/F11
1695  print 1004,tb,f11,f22,f33,f44,f12,f34
1696  f(i1,1,1)=f11
1697  f(i1,1,2)=f12
1698  f(i1,1,3)=0
1699  f(i1,1,4)=0
1700  f(i1,2,1)=f12
1701  f(i1,2,2)=f22
1702  f(i1,2,3)=0
1703  f(i1,2,4)=0
1704  f(i1,3,1)=0
1705  f(i1,3,2)=0
1706  f(i1,3,3)=f33
1707  f(i1,3,4)=f34
1708  f(i1,4,1)=0
1709  f(i1,4,2)=0
1710  f(i1,4,3)=-f34
1711  f(i1,4,4)=f44
1712  500 CONTINUE
1713  print 1000
1714  1004 FORMAT(' ',f6.2,6f11.4)
1715  RETURN
1716  END
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
double precision function zeroin(AX, BX, F, TOL)
Definition: tmd.lp.f:1395
subroutine hovenr(L1, A1, A2, A3, A4, B1, B2)
Definition: tmd.lp.f:1561
subroutine ccg(N, N1, NMAX, K1, K2, GG)
Definition: tmd.lp.f:1209
subroutine tmatr0(NGAUSS, X, W, AN, ANN, S, SS, PPI, PIR, PII, R, DR, DDR, DRR, DRI, NMAX, NCHECK)
Definition: ampld.lp.f:1301
subroutine gsp(NMAX, CSCA, LAM, ALF1, ALF2, ALF3, ALF4, BET1, BET2, LMAX)
Definition: tmd.lp.f:813
float f3(float z)
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
subroutine const(NGAUSS, NMAX, MMAX, P, X, W, AN, ANN, S, SS, NP, EPS)
Definition: ampld.lp.f:924
subroutine signum
Definition: tmd.lp.f:1178
#define fac
float f2(float y)
subroutine matr(A1, A2, A3, A4, B1, B2, LMAX, NPNA, F)
Definition: tmd.lp.f:1615
subroutine sareac(EPS, RAT)
Definition: ampld.lp.f:1951
subroutine sarea(D, RAT)
Definition: ampld.lp.f:1902
subroutine surfch(N, E, RAT)
Definition: ampld.lp.f:1920
#define max(A, B)
Definition: main_biosmap.c:61
===========================================================================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
#define real
subroutine fact
Definition: tmd.lp.f:1161
#define min(A, B)
Definition: main_biosmap.c:62
subroutine power(A, B, R1, R2)
Definition: tmd.lp.f:1380
subroutine ccgin(N, N1, M, MM, G)
Definition: tmd.lp.f:1318
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
subroutine distrb(NNK, YY, WY, NDISTR, AA, BB, GAM, R1, R2, REFF, VEFF, PI)
Definition: tmd.lp.f:1474
#define f
Definition: l1_czcs_hdf.c:701
Definition: L3File.cpp:10
subroutine gauss(x1, x2, x, w, n)
Definition: 6sm1.f:4055
#define abs(a)
Definition: misc.h:90
subroutine direct(N, M, N1, M1, NN, MM, C)
Definition: tmd.lp.f:1299