OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
tm_ampld.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 THE AMPLITUDE AND PHASE MATRICES FOR
14 C A PARTICLE WITH AN AXIALLY SYMMETRIC SHAPE
15 
16 C This version of the code uses DOUBLE PRECISION variables,
17 C is applicable to spheroids, finite circular cylinders,
18 C Chebyshev particles, and generalized Chebyshev particles
19 C (distorted water drops), and must be used along with the
20 C accompanying files lpd.f and tm_ampld.par.f.
21 
22 C Last update 08/06/2005
23 
24 C The code has been developed by Michael Mishchenko at the NASA
25 C Goddard Institute for Space Studies, New York. The development
26 C of the code was supported by the NASA Radiation Sciences Program.
27 
28 C We only request that in any publication using the code the source of
29 C the code be acknowledged and relevant references (see below) be
30 C made.
31 
32 C The computational method is based on the T-matrix approach
33 C [P. C. Waterman, Phys. Rev. D 3, 825 (1971)], also known as
34 C the extended boundary condition method. The method was
35 C improved in the following papers:
36 
37 C 1. M. I. Mishchenko and L. D. Travis, T-matrix computations
38 C of light scattering by large spheroidal particles,
39 C Opt. Commun., vol. 109, 16-21 (1994).
40 C
41 C 2. M. I. Mishchenko, L. D. Travis, and A. Macke, Scattering
42 C of light by polydisperse, randomly oriented, finite
43 C circular cylinders, Appl. Opt., vol. 35, 4927-4940 (1996).
44 C
45 C 3. D. J. Wielaard, M. I. Mishchenko, A. Macke, and B. E. Carlson,
46 C Improved T-matrix computations for large, nonabsorbing and
47 C weakly absorbing nonspherical particles and comparison
48 C with geometrical optics approximation, Appl. Opt., vol. 36,
49 C 4305-4313 (1997).
50 C
51 C A general review of the T-matrix approach can be found in
52 C
53 C 4. M. I. Mishchenko, L. D. Travis, and D. W. Mackowski,
54 C T-matrix computations of light scattering by nonspherical
55 C particles: a review, J. Quant. Spectrosc. Radiat.
56 C Transfer, vol. 55, 535-575 (1996).
57 C
58 C Additional useful information is contained in the paper
59 C
60 C 5. M. I. Mishchenko and L. D. Travis, Capabilities and
61 C limitations of a current FORTRAN implementation of the
62 C T-matrix method for randomly oriented, rotationally
63 C symmetric scatterers, J. Quant. Spectrosc. Radiat. Transfer,
64 C vol. 60, 309-324 (1998).
65 C
66 C The definitions and notation used can also be found in
67 C
68 C 6. M. I. Mishchenko, Calculation of the amplitude matrix
69 C for a nonspherical particle in a fixed orientation,
70 C Appl. Opt. vol. 39, 1026-1031 (2000).
71 
72 C These papers are available in the .pdf format at the web site
73 C
74 C http://www.giss.nasa.gov/~crmim/publications/
75 C
76 C or in hardcopy upon request from Michael Mishchenko
77 C Please e-mail your request to crmim@giss.nasa.gov.
78 C
79 C A comprehensive book "Scattering, Absorption, and Emission of
80 C Light by Small Particles" (Cambridge University Press, Cambridge,
81 C 2002) is also available in the .pdf format at the web site
82 C
83 C http://www.giss.nasa.gov/~crmim/books.html
84 
85 C One of the main advantages of the T-matrix method is that the
86 C T-matrix for a given nonspherical particle needs to be computed
87 C only once and then can be used any number of times for computing
88 C the amplitude and phase matrices for any directions of the incident
89 C and scattered beams. This makes the T-matrix method
90 C extremely convenient and efficient in averaging over particle
91 C orientations and/or directions of incidence (or scattering).
92 
93 C The use of extended precision variables (Ref. 1) can significantly
94 C increase the maximum convergent equivalent-sphere size parameter
95 C and make it well larger than 100 (depending on refractive index
96 C and aspect ratio). The extended-precision code is also available.
97 C However, the use of extended precision varibales results in a
98 C greater consumption of CPU time.
99 C On IBM RISC workstations, that code is approximately
100 C five times slower than this double-precision code. The
101 C CPU time difference between the double-precision and extended-
102 C precision codes can be larger on supercomputers.
103 
104 C This is the first part of the full T-matrix code. The second part,
105 C lpq.f, is completely independent of the first part. It contains no
106 C T-matrix-specific subroutines and can be compiled separately.
107 C The second part of the code replaces the previously implemented
108 C standard matrix inversion scheme based on Gaussian elimination
109 C by a scheme based on the LU factorization technique.
110 C As described in Ref. 3 above, the use of the LU factorization is
111 C especially beneficial for nonabsorbing or weakly absorbing particles.
112 C In this code we use the LAPACK implementation of the LU factorization
113 C scheme. LAPACK stands for Linear Algebra PACKage. The latter is
114 C publicly available at the following internet site:
115 C
116 C http://www.netlib.org/lapack/
117 
118 
119 C INPUT PARAMETERS:
120 C
121 C AXI - equivalent-sphere radius
122 C RAT = 1 - particle size is specified in terms of the
123 C equal-volume-sphere radius
124 C RAT.ne.1 - particle size is specified in terms of the
125 C equal-surface-area-sphere radius
126 C LAM - WAVELENGTH OF INCIDENT LIGHT
127 C MRR and MRI - real and imaginary parts of the refractive
128 C index
129 C EPS and NP - specify the shape of the particles.
130 C For spheroids NP=-1 and EPS is the ratio of the
131 C horizontal to rotational axes. EPS is larger than
132 C 1 for oblate spheroids and smaller than 1 for
133 C prolate spheroids.
134 C For cylinders NP=-2 and EPS is the ratio of the
135 C diameter to the length.
136 C For Chebyshev particles NP must be positive and
137 C is the degree of the Chebyshev polynomial, while
138 C EPS is the deformation parameter (Ref. 5).
139 C For generalized Chebyshev particles (describing the shape
140 C of distorted water drops) NP=-3. The coefficients
141 C of the Chebyshev polynomial expansion of the particle
142 C shape (Ref. 7) are specified in subroutine DROP.
143 C DDELT - accuracy of the computations
144 C NDGS - parameter controlling the number of division points
145 C in computing integrals over the particle surface (Ref. 5).
146 C For compact particles, the recommended value is 2.
147 C For highly aspherical particles larger values (3, 4,...)
148 C may be necessary to obtain convergence.
149 C The code does not check convergence over this parameter.
150 C Therefore, control comparisons of results obtained with
151 C different NDGS-values are recommended.
152 C ALPHA and BETA - Euler angles (in degrees) specifying the orientation
153 C of the scattering particle relative to the laboratory reference
154 C frame (Refs. 6 and 7).
155 C THET0 - zenith angle of the incident beam in degrees
156 C THET - zenith angle of the scattered beam in degrees
157 C PHI0 - azimuth angle of the incident beam in degrees
158 C PHI - azimuth angle of the scattered beam in degrees
159 C (Refs. 6 and 7)
160 C ALPHA, BETA, THET0, THET, PHI0, and PHI are specified at the end of
161 C the main program before the line
162 C
163 C "CALL AMPL (NMAX,...)"
164 C
165 C The part of the main program following the line
166 C
167 C "COMPUTATION OF THE AMPLITUDE AND PHASE MATRICES"
168 C
169 C can be repeated any number of times. At this point the T-matrix
170 C for the given scattering particle has already
171 C been fully computed and can be repeatedly used in computations
172 C for any directions of illumination and scattering and any particle
173 C orientations.
174 
175 C OUTPUT PARAMETERS:
176 C
177 C Elements of the 2x2 amplitude matrix
178 C Elements of the 4x4 phase matrix
179 
180 C Note that LAM and AXI must be given in the same units of length
181 C (e.g., microns).
182 
183 C The convergence of the T-matrix method for particles with
184 C different sizes, refractive indices, and aspect ratios can be
185 C dramatically different. Usually, large sizes and large aspect
186 C ratios cause problems. The user of this code
187 C should "play" a little with different input parameters in
188 C order to get an idea of the range of applicability of this
189 C technique. Sometimes decreasing the aspect ratio
190 C from 3 to 2 can increase the maximum convergent equivalent-
191 C sphere size parameter by a factor of several.
192 
193 C The CPU time required rapidly increases with increasing ratio
194 C radius/wavelength and/or with increasing particle asphericity.
195 C This should be taken into account in planning massive computations.
196 
197 C Execution can be automatically terminated if dimensions of certain
198 C arrays are not big enough or if the convergence procedure decides
199 C that the accuracy of double-precision variables is insufficient
200 C to obtain a converged T-matrix solution for given particles.
201 C In all cases, a message appears explaining
202 C the cause of termination.
203 
204 C The message
205 C "WARNING: W IS GREATER THAN 1"
206 C means that the single-scattering albedo exceeds the maximum
207 C possible value 1. If W is greater than 1 by more than
208 C DDELT, this message can be an indication of numerical
209 C instability caused by extreme values of particle parameters.
210 
211 C The message "WARNING: NGAUSS=NPNG1" means that convergence over
212 C the parameter NG (see Ref. 2) cannot be obtained for the NPNG1
213 C value specified in the PARAMETER statement in the file tm_ampld.par.f.
214 C Often this is not a serious problem, especially for compact
215 C particles.
216 
217 C Larger and/or more aspherical particles may require larger
218 C values of the parameters NPN1, NPN4, and NPNG1 in the file
219 C tm_ampld.par.f. It is recommended to keep NPN1=NPN4+25 and
220 C NPNG1=3*NPN1. Note that the memory requirement increases
221 C as the third power of NPN4. If the memory of
222 C a computer is too small to accomodate the code in its current
223 C setting, the parameters NPN1, NPN4, and NPNG1 should be
224 C decreased. However, this will decrease the maximum size parameter
225 C that can be handled by the code.
226 
227 C In some cases any increases of NPN1 will not make the T-matrix
228 C computations convergent. This means that the particle is just
229 C too "bad" (extreme size parameter and/or extreme aspect ratio
230 C and/or extreme refractive index).
231 C The main program contains several PRINT statements which are
232 C currently commentd out. If uncommented, these statements will
233 C produce numbers which show the convergence rate and can be
234 C used to determine whether T-matrix computations for given particle
235 C parameters will converge at all, whatever the parameter NPN1 is.
236 
237 C Some of the common blocks are used to save memory rather than
238 C to transfer data. Therefore, if a compiler produces a warning
239 C message that the lengths of a common block are different in
240 C different subroutines, this is not a real problem.
241 
242 C The recommended value of DDELT is 0.001. For bigger values,
243 C false convergence can be obtained.
244 
245 C In computations for spheres, use EPS=1.000001 instead of EPS=1.
246 C EPS=1 can cause overflows in some rare cases.
247 
248 C For some compilers, DACOS must be raplaced by DARCOS and DASIN
249 C by DARSIN.
250 
251 C I would highly appreciate informing me of any problems encountered
252 C with this code. Please send your message to the following
253 C e-mail address: CRMIM@GISS.NASA.GOV.
254 
255 C WHILE THE COMPUTER PROGRAM HAS BEEN TESTED FOR A VARIETY OF CASES,
256 C IT IS NOT INCONCEIVABLE THAT IT CONTAINS UNDETECTED ERRORS. ALSO,
257 C INPUT PARAMETERS CAN BE USED WHICH ARE OUTSIDE THE ENVELOPE OF
258 C VALUES FOR WHICH RESULTS ARE COMPUTED ACCURATELY. FOR THIS REASON,
259 C THE AUTHORS AND THEIR ORGANIZATION DISCLAIM ALL LIABILITY FOR
260 C ANY DAMAGES THAT MAY RESULT FROM THE USE OF THE PROGRAM.
261 
262  SUBROUTINE calctmat(AXI,RAT,LAM,MRR,MRI,EPS,NP,DDELT,NDGS,NMAX)
263  IMPLICIT REAL*8 (a-h,o-z)
264 
265 Cf2py intent(in) axi
266 Cf2py intent(in) rat
267 Cf2py intent(in) lam
268 Cf2py intent(in) mrr
269 Cf2py intent(in) mri
270 Cf2py intent(in) eps
271 Cf2py intent(in) np
272 Cf2py intent(in) ddelt
273 Cf2py intent(in) ndgs
274 Cf2py intent(out) nmax
275 
276  include 'tm_tmd.par.f'
277  real*8 lam,mrr,mri,x(npng2),w(npng2),s(npng2),ss(npng2),
278  * an(npn1),r(npng2),dr(npng2),
279  * ddr(npng2),drr(npng2),dri(npng2),ann(npn1,npn1)
280  real*8 tr1(npn2,npn2),ti1(npn2,npn2)
281  real*8 xalpha(300),xbeta(300),walpha(300),wbeta(300)
282  real*4
283  & rt11(npn6,npn4,npn4),rt12(npn6,npn4,npn4),
284  & rt21(npn6,npn4,npn4),rt22(npn6,npn4,npn4),
285  & it11(npn6,npn4,npn4),it12(npn6,npn4,npn4),
286  & it21(npn6,npn4,npn4),it22(npn6,npn4,npn4)
287 
288  COMMON /ct/ tr1,ti1
289  COMMON /tmatc/ rt11,rt12,rt21,rt22,it11,it12,it21,it22
290  COMMON /choice/ ichoice
291 
292 C OPEN FILES *******************************************************
293 
294 C OPEN (6,FILE='test')
295 
296 C INPUT DATA ********************************************************
297 
298 C AXI=10D0
299 C RAT=0.1 D0
300 C LAM=DACOS(-1D0)*2D0
301 C MRR=1.5 D0
302 C MRI=0.02 D0
303 C EPS=0.5 D0
304 C NP=-1
305 C DDELT=0.001D0
306 C NDGS=2
307 
308  p=dacos(-1d0)
309  ncheck=0
310  IF (np.EQ.-1.OR.np.EQ.-2) ncheck=1
311  IF (np.GT.0.AND.(-1)**np.EQ.1) ncheck=1
312 C WRITE (6,5454) NCHECK
313 C 5454 FORMAT ('NCHECK=',I1)
314  IF (dabs(rat-1d0).GT.1d-8.AND.np.EQ.-1) CALL sarea (eps,rat)
315  IF (dabs(rat-1d0).GT.1d-8.AND.np.GE.0) CALL surfch(np,eps,rat)
316  IF (dabs(rat-1d0).GT.1d-8.AND.np.EQ.-2) CALL sareac (eps,rat)
317  IF (np.EQ.-3) CALL drop (rat)
318 C PRINT 8000, RAT
319 C 8000 FORMAT ('RAT=',F8.6)
320 C IF(NP.EQ.-1.AND.EPS.GE.1D0) PRINT 7000,EPS
321 C IF(NP.EQ.-1.AND.EPS.LT.1D0) PRINT 7001,EPS
322 C IF(NP.GE.0) PRINT 7100,NP,EPS
323 C IF(NP.EQ.-2.AND.EPS.GE.1D0) PRINT 7150,EPS
324 C IF(NP.EQ.-2.AND.EPS.LT.1D0) PRINT 7151,EPS
325 C IF(NP.EQ.-3) PRINT 7160
326 C PRINT 7400, LAM,MRR,MRI
327 C PRINT 7200,DDELT
328 C 7000 FORMAT('OBLATE SPHEROIDS, A/B=',F11.7)
329 C 7001 FORMAT('PROLATE SPHEROIDS, A/B=',F11.7)
330 C 7100 FORMAT('CHEBYSHEV PARTICLES, T',
331 C & I1,'(',F5.2,')')
332 C 7150 FORMAT('OBLATE CYLINDERS, D/L=',F11.7)
333 C 7151 FORMAT('PROLATE CYLINDERS, D/L=',F11.7)
334 C 7160 FORMAT('GENERALIZED CHEBYSHEV PARTICLES')
335 C 7200 FORMAT ('ACCURACY OF COMPUTATIONS DDELT = ',D8.2)
336 C 7400 FORMAT('LAM=',F10.6,3X,'MRR=',D10.4,3X,'MRI=',D10.4)
337 C DDELT=0.1D0*DDELT
338 C IF (DABS(RAT-1D0).LE.1D-6) PRINT 8003, AXI
339 C IF (DABS(RAT-1D0).GT.1D-6) PRINT 8004, AXI
340 C 8003 FORMAT('EQUAL-VOLUME-SPHERE RADIUS=',F8.4)
341 C 8004 FORMAT('EQUAL-SURFACE-AREA-SPHERE RADIUS=',F8.4)
342  a=rat*axi
343  xev=2d0*p*a/lam
344  ixxx=xev+4.05d0*xev**0.333333d0
345  inm1=max0(4,ixxx)
346 C IF (INM1.GE.NPN1) PRINT 7333, NPN1
347  IF (inm1.GE.npn1) stop
348 C 7333 FORMAT('CONVERGENCE IS NOT OBTAINED FOR NPN1=',I3,
349 C & '. EXECUTION TERMINATED')
350  qext1=0d0
351  qsca1=0d0
352  DO 50 nma=inm1,npn1
353  nmax=nma
354  mmax=1
355  ngauss=nmax*ndgs
356 C IF (NGAUSS.GT.NPNG1) PRINT 7340, NGAUSS
357  IF (ngauss.GT.npng1) stop
358 C 7340 FORMAT('NGAUSS =',I3,' I.E. IS GREATER THAN NPNG1.',
359 C & ' EXECUTION TERMINATED')
360 C 7334 FORMAT(' NMAX =', I3,' DC2=',D8.2,' DC1=',D8.2)
361  CALL const(ngauss,nmax,mmax,p,x,w,an,ann,s,ss,np,eps)
362  CALL vary(lam,mrr,mri,a,eps,np,ngauss,x,p,ppi,pir,pii,r,
363  & dr,ddr,drr,dri,nmax)
364  CALL tmatr0 (ngauss,x,w,an,ann,s,ss,ppi,pir,pii,r,dr,
365  & ddr,drr,dri,nmax,ncheck)
366  qext=0d0
367  qsca=0d0
368  DO 4 n=1,nmax
369  n1=n+nmax
370  tr1nn=tr1(n,n)
371  ti1nn=ti1(n,n)
372  tr1nn1=tr1(n1,n1)
373  ti1nn1=ti1(n1,n1)
374  dn1=dfloat(2*n+1)
375  qsca=qsca+dn1*(tr1nn*tr1nn+ti1nn*ti1nn
376  & +tr1nn1*tr1nn1+ti1nn1*ti1nn1)
377  qext=qext+(tr1nn+tr1nn1)*dn1
378  4 CONTINUE
379  dsca=dabs((qsca1-qsca)/qsca)
380  dext=dabs((qext1-qext)/qext)
381  qext1=qext
382  qsca1=qsca
383 C PRINT 7334, NMAX,DSCA,DEXT
384  IF(dsca.LE.ddelt.AND.dext.LE.ddelt) GO TO 55
385 C IF (NMA.EQ.NPN1) PRINT 7333, NPN1
386  IF (nma.EQ.npn1) stop
387  50 CONTINUE
388  55 nnnggg=ngauss+1
389  mmax=nmax
390 C IF (NGAUSS.EQ.NPNG1) PRINT 7336
391  IF (ngauss.EQ.npng1) GO TO 155
392  DO 150 ngaus=nnnggg,npng1
393  ngauss=ngaus
394  nggg=2*ngauss
395 C 7336 FORMAT('WARNING: NGAUSS=NPNG1')
396 C 7337 FORMAT(' NG=',I3,' DC2=',D8.2,' DC1=',D8.2)
397  CALL const(ngauss,nmax,mmax,p,x,w,an,ann,s,ss,np,eps)
398  CALL vary(lam,mrr,mri,a,eps,np,ngauss,x,p,ppi,pir,pii,r,
399  & dr,ddr,drr,dri,nmax)
400  CALL tmatr0 (ngauss,x,w,an,ann,s,ss,ppi,pir,pii,r,dr,
401  & ddr,drr,dri,nmax,ncheck)
402  qext=0d0
403  qsca=0d0
404  DO 104 n=1,nmax
405  n1=n+nmax
406  tr1nn=tr1(n,n)
407  ti1nn=ti1(n,n)
408  tr1nn1=tr1(n1,n1)
409  ti1nn1=ti1(n1,n1)
410  dn1=dfloat(2*n+1)
411  qsca=qsca+dn1*(tr1nn*tr1nn+ti1nn*ti1nn
412  & +tr1nn1*tr1nn1+ti1nn1*ti1nn1)
413  qext=qext+(tr1nn+tr1nn1)*dn1
414  104 CONTINUE
415  dsca=dabs((qsca1-qsca)/qsca)
416  dext=dabs((qext1-qext)/qext)
417 C PRINT 7337, NGGG,DSCA,DEXT
418  qext1=qext
419  qsca1=qsca
420  IF(dsca.LE.ddelt.AND.dext.LE.ddelt) GO TO 155
421 C IF (NGAUS.EQ.NPNG1) PRINT 7336
422  150 CONTINUE
423  155 CONTINUE
424  qsca=0d0
425  qext=0d0
426  nnm=nmax*2
427  DO 204 n=1,nnm
428  qext=qext+tr1(n,n)
429  204 CONTINUE
430  DO 213 n2=1,nmax
431  nn2=n2+nmax
432  DO 213 n1=1,nmax
433  nn1=n1+nmax
434  zz1=tr1(n1,n2)
435  rt11(1,n1,n2)=zz1
436  zz2=ti1(n1,n2)
437  it11(1,n1,n2)=zz2
438  zz3=tr1(n1,nn2)
439  rt12(1,n1,n2)=zz3
440  zz4=ti1(n1,nn2)
441  it12(1,n1,n2)=zz4
442  zz5=tr1(nn1,n2)
443  rt21(1,n1,n2)=zz5
444  zz6=ti1(nn1,n2)
445  it21(1,n1,n2)=zz6
446  zz7=tr1(nn1,nn2)
447  rt22(1,n1,n2)=zz7
448  zz8=ti1(nn1,nn2)
449  it22(1,n1,n2)=zz8
450  qsca=qsca+zz1*zz1+zz2*zz2+zz3*zz3+zz4*zz4
451  & +zz5*zz5+zz6*zz6+zz7*zz7+zz8*zz8
452  213 CONTINUE
453 C PRINT 7800,0,DABS(QEXT),QSCA,NMAX
454  DO 220 m=1,nmax
455  CALL tmatr(m,ngauss,x,w,an,ann,s,ss,ppi,pir,pii,r,dr,
456  & ddr,drr,dri,nmax,ncheck)
457  nm=nmax-m+1
458  m1=m+1
459  qsc=0d0
460  DO 214 n2=1,nm
461  nn2=n2+m-1
462  n22=n2+nm
463  DO 214 n1=1,nm
464  nn1=n1+m-1
465  n11=n1+nm
466  zz1=tr1(n1,n2)
467  rt11(m1,nn1,nn2)=zz1
468  zz2=ti1(n1,n2)
469  it11(m1,nn1,nn2)=zz2
470  zz3=tr1(n1,n22)
471  rt12(m1,nn1,nn2)=zz3
472  zz4=ti1(n1,n22)
473  it12(m1,nn1,nn2)=zz4
474  zz5=tr1(n11,n2)
475  rt21(m1,nn1,nn2)=zz5
476  zz6=ti1(n11,n2)
477  it21(m1,nn1,nn2)=zz6
478  zz7=tr1(n11,n22)
479  rt22(m1,nn1,nn2)=zz7
480  zz8=ti1(n11,n22)
481  it22(m1,nn1,nn2)=zz8
482  qsc=qsc+(zz1*zz1+zz2*zz2+zz3*zz3+zz4*zz4
483  & +zz5*zz5+zz6*zz6+zz7*zz7+zz8*zz8)*2d0
484  214 CONTINUE
485  nnm=2*nm
486  qxt=0d0
487  DO 215 n=1,nnm
488  qxt=qxt+tr1(n,n)*2d0
489  215 CONTINUE
490  qsca=qsca+qsc
491  qext=qext+qxt
492 C PRINT 7800,M,DABS(QXT),QSC,NMAX
493 C 7800 FORMAT(' m=',I3,' qxt=',D12.6,' qsc=',D12.6,
494 C & ' nmax=',I3)
495  220 CONTINUE
496  walb=-qsca/qext
497 C IF (WALB.GT.1D0+DDELT) PRINT 9111
498 C 9111 FORMAT ('WARNING: W IS GREATER THAN 1')
499  RETURN
500  END
501 
502 C COMPUTATION OF THE AMPLITUDE AND PHASE MATRICES
503 
504 C ALPHA=145D0
505 C BETA=52D0
506 C THET0=56D0
507 C THET=65D0
508 C PHI0=114D0
509 C PHI=128D0
510 
511  SUBROUTINE calcampl(NMAX,LAM,THET0,THET,PHI0,PHI,ALPHA,BETA,S,Z)
512  IMPLICIT REAL*8 (a-h,o-z)
513  real*8 lam
514  COMPLEX*16 S(2,2)
515  real*8 z(4,4)
516 Cf2py intent(in) nmax
517 Cf2py intent(in) lam
518 Cf2py intent(in) thet0
519 Cf2py intent(in) thet
520 Cf2py intent(in) phi0
521 Cf2py intent(in) phi
522 Cf2py intent(in) alpha
523 Cf2py intent(in) beta
524 Cf2py intent(out) s
525 Cf2py intent(out) z
526 
527 C AMPLITUDE MATRIX [Eqs. (2)-(4) of Ref. 6]
528  CALL ampl (nmax,lam,thet0,thet,phi0,phi,alpha,beta,
529  & s(1,1),s(1,2),s(2,1),s(2,2))
530 C PHASE MATRIX [Eqs. (13)-(29) of Ref. 6]
531  z(1,1)=0.5d0*(s(1,1)*dconjg(s(1,1))+s(1,2)*dconjg(s(1,2))
532  & +s(2,1)*dconjg(s(2,1))+s(2,2)*dconjg(s(2,2)))
533  z(1,2)=0.5d0*(s(1,1)*dconjg(s(1,1))-s(1,2)*dconjg(s(1,2))
534  & +s(2,1)*dconjg(s(2,1))-s(2,2)*dconjg(s(2,2)))
535  z(1,3)=-s(1,1)*dconjg(s(1,2))-s(2,2)*dconjg(s(2,1))
536  z(1,4)=(0d0,1d0)*(s(1,1)*dconjg(s(1,2))-s(2,2)*dconjg(s(2,1)))
537  z(2,1)=0.5d0*(s(1,1)*dconjg(s(1,1))+s(1,2)*dconjg(s(1,2))
538  & -s(2,1)*dconjg(s(2,1))-s(2,2)*dconjg(s(2,2)))
539  z(2,2)=0.5d0*(s(1,1)*dconjg(s(1,1))-s(1,2)*dconjg(s(1,2))
540  & -s(2,1)*dconjg(s(2,1))+s(2,2)*dconjg(s(2,2)))
541  z(2,3)=-s(1,1)*dconjg(s(1,2))+s(2,2)*dconjg(s(2,1))
542  z(2,4)=(0d0,1d0)*(s(1,1)*dconjg(s(1,2))+s(2,2)*dconjg(s(2,1)))
543  z(3,1)=-s(1,1)*dconjg(s(2,1))-s(2,2)*dconjg(s(1,2))
544  z(3,2)=-s(1,1)*dconjg(s(2,1))+s(2,2)*dconjg(s(1,2))
545  z(3,3)=s(1,1)*dconjg(s(2,2))+s(1,2)*dconjg(s(2,1))
546  z(3,4)=(0d0,-1d0)*(s(1,1)*dconjg(s(2,2))+s(2,1)*dconjg(s(1,2)))
547  z(4,1)=(0d0,1d0)*(s(2,1)*dconjg(s(1,1))+s(2,2)*dconjg(s(1,2)))
548  z(4,2)=(0d0,1d0)*(s(2,1)*dconjg(s(1,1))-s(2,2)*dconjg(s(1,2)))
549  z(4,3)=(0d0,-1d0)*(s(2,2)*dconjg(s(1,1))-s(1,2)*dconjg(s(2,1)))
550  z(4,4)=s(2,2)*dconjg(s(1,1))-s(1,2)*dconjg(s(2,1))
551 C WRITE (6,5000)
552 C WRITE (6,5001) Z11,Z12,Z13,Z14
553 C WRITE (6,5001) Z21,Z22,Z23,Z24
554 C WRITE (6,5001) Z31,Z32,Z33,Z34
555 C WRITE (6,5001) Z41,Z42,Z43,Z44
556 
557 C 5000 FORMAT ('PHASE MATRIX')
558 C 5001 FORMAT (4F10.4)
559 
560 C ITIME=MCLOCK()
561 C TIME=DFLOAT(ITIME)/6000D0
562 C PRINT 1001,TIME
563 C 1001 FORMAT (' time =',F8.2,' min')
564 C STOP
565  RETURN
566  END
567 
568 C********************************************************************
569 
570 C CALCULATION OF THE AMPLITUDE MATRIX
571 
572  SUBROUTINE ampl (NMAX,DLAM,TL,TL1,PL,PL1,ALPHA,BETA,
573  & VV,VH,HV,HH)
574  include 'tm_tmd.par.f'
575  IMPLICIT REAL*8 (a-b,d-h,o-z), COMPLEX*16 (C)
576  real*8 al(3,2),al1(3,2),ap(2,3),ap1(2,3),b(3,3),
577  * r(2,2),r1(2,2),c(3,2),ca,cb,ct,cp,ctp,cpp,ct1,cp1,
578  * ctp1,cpp1
579  real*8 dv1(npn6),dv2(npn6),dv01(npn6),dv02(npn6)
580  real*4
581  & tr11(npn6,npn4,npn4),tr12(npn6,npn4,npn4),
582  & tr21(npn6,npn4,npn4),tr22(npn6,npn4,npn4),
583  & ti11(npn6,npn4,npn4),ti12(npn6,npn4,npn4),
584  & ti21(npn6,npn4,npn4),ti22(npn6,npn4,npn4)
585  COMPLEX*16 CAL(NPN4,NPN4),VV,VH,HV,HH
586  COMMON /tmatc/ tr11,tr12,tr21,tr22,ti11,ti12,ti21,ti22
587 
588  IF (alpha.LT.0d0.OR.alpha.GT.360d0.OR.
589  & beta.LT.0d0.OR.beta.GT.180d0.OR.
590  & tl.LT.0d0.OR.tl.GT.180d0.OR.
591  & tl1.LT.0d0.OR.tl1.GT.180d0.OR.
592  & pl.LT.0d0.OR.pl.GT.360d0.OR.
593  & pl1.LT.0d0.OR.pl1.GT.360d0) THEN
594 C WRITE (6,2000)
595  stop
596  ELSE
597  CONTINUE
598  ENDIF
599 C 2000 FORMAT ('AN ANGULAR PARAMETER IS OUTSIDE ITS',
600 C & ' ALLOWABLE RANGE')
601  pin=dacos(-1d0)
602  pin2=pin*0.5d0
603  pi=pin/180d0
604  alph=alpha*pi
605  bet=beta*pi
606  thetl=tl*pi
607  phil=pl*pi
608  thetl1=tl1*pi
609  phil1=pl1*pi
610 
611  eps=1d-7
612  IF (thetl.LT.pin2) thetl=thetl+eps
613  IF (thetl.GT.pin2) thetl=thetl-eps
614  IF (thetl1.LT.pin2) thetl1=thetl1+eps
615  IF (thetl1.GT.pin2) thetl1=thetl1-eps
616  IF (phil.LT.pin) phil=phil+eps
617  IF (phil.GT.pin) phil=phil-eps
618  IF (phil1.LT.pin) phil1=phil1+eps
619  IF (phil1.GT.pin) phil1=phil1-eps
620  IF (bet.LE.pin2.AND.pin2-bet.LE.eps) bet=bet-eps
621  IF (bet.GT.pin2.AND.bet-pin2.LE.eps) bet=bet+eps
622 
623 C_____________COMPUTE THETP, PHIP, THETP1, AND PHIP1, EQS. (8), (19), AND (20)
624 
625  cb=dcos(bet)
626  sb=dsin(bet)
627  ct=dcos(thetl)
628  st=dsin(thetl)
629  cp=dcos(phil-alph)
630  sp=dsin(phil-alph)
631  ctp=ct*cb+st*sb*cp
632  thetp=dacos(ctp)
633  cpp=cb*st*cp-sb*ct
634  spp=st*sp
635  phip=datan(spp/cpp)
636  IF (phip.GT.0d0.AND.sp.LT.0d0) phip=phip+pin
637  IF (phip.LT.0d0.AND.sp.GT.0d0) phip=phip+pin
638  IF (phip.LT.0d0) phip=phip+2d0*pin
639 
640  ct1=dcos(thetl1)
641  st1=dsin(thetl1)
642  cp1=dcos(phil1-alph)
643  sp1=dsin(phil1-alph)
644  ctp1=ct1*cb+st1*sb*cp1
645  thetp1=dacos(ctp1)
646  cpp1=cb*st1*cp1-sb*ct1
647  spp1=st1*sp1
648  phip1=datan(spp1/cpp1)
649  IF (phip1.GT.0d0.AND.sp1.LT.0d0) phip1=phip1+pin
650  IF (phip1.LT.0d0.AND.sp1.GT.0d0) phip1=phip1+pin
651  IF (phip1.LT.0d0) phip1=phip1+2d0*pin
652 
653 C____________COMPUTE MATRIX BETA, EQ. (21)
654 
655  ca=dcos(alph)
656  sa=dsin(alph)
657  b(1,1)=ca*cb
658  b(1,2)=sa*cb
659  b(1,3)=-sb
660  b(2,1)=-sa
661  b(2,2)=ca
662  b(2,3)=0d0
663  b(3,1)=ca*sb
664  b(3,2)=sa*sb
665  b(3,3)=cb
666 
667 C____________COMPUTE MATRICES AL AND AL1, EQ. (14)
668 
669  cp=dcos(phil)
670  sp=dsin(phil)
671  cp1=dcos(phil1)
672  sp1=dsin(phil1)
673  al(1,1)=ct*cp
674  al(1,2)=-sp
675  al(2,1)=ct*sp
676  al(2,2)=cp
677  al(3,1)=-st
678  al(3,2)=0d0
679  al1(1,1)=ct1*cp1
680  al1(1,2)=-sp1
681  al1(2,1)=ct1*sp1
682  al1(2,2)=cp1
683  al1(3,1)=-st1
684  al1(3,2)=0d0
685 
686 C____________COMPUTE MATRICES AP^(-1) AND AP1^(-1), EQ. (15)
687 
688  ct=ctp
689  st=dsin(thetp)
690  cp=dcos(phip)
691  sp=dsin(phip)
692  ct1=ctp1
693  st1=dsin(thetp1)
694  cp1=dcos(phip1)
695  sp1=dsin(phip1)
696  ap(1,1)=ct*cp
697  ap(1,2)=ct*sp
698  ap(1,3)=-st
699  ap(2,1)=-sp
700  ap(2,2)=cp
701  ap(2,3)=0d0
702  ap1(1,1)=ct1*cp1
703  ap1(1,2)=ct1*sp1
704  ap1(1,3)=-st1
705  ap1(2,1)=-sp1
706  ap1(2,2)=cp1
707  ap1(2,3)=0d0
708 
709 C____________COMPUTE MATRICES R AND R^(-1), EQ. (13)
710  DO i=1,3
711  DO j=1,2
712  x=0d0
713  DO k=1,3
714  x=x+b(i,k)*al(k,j)
715  ENDDO
716  c(i,j)=x
717  ENDDO
718  ENDDO
719  DO i=1,2
720  DO j=1,2
721  x=0d0
722  DO k=1,3
723  x=x+ap(i,k)*c(k,j)
724  ENDDO
725  r(i,j)=x
726  ENDDO
727  ENDDO
728  DO i=1,3
729  DO j=1,2
730  x=0d0
731  DO k=1,3
732  x=x+b(i,k)*al1(k,j)
733  ENDDO
734  c(i,j)=x
735  ENDDO
736  ENDDO
737  DO i=1,2
738  DO j=1,2
739  x=0d0
740  DO k=1,3
741  x=x+ap1(i,k)*c(k,j)
742  ENDDO
743  r1(i,j)=x
744  ENDDO
745  ENDDO
746  d=1d0/(r1(1,1)*r1(2,2)-r1(1,2)*r1(2,1))
747  x=r1(1,1)
748  r1(1,1)=r1(2,2)*d
749  r1(1,2)=-r1(1,2)*d
750  r1(2,1)=-r1(2,1)*d
751  r1(2,2)=x*d
752 
753  ci=(0d0,1d0)
754  DO 5 nn=1,nmax
755  DO 5 n=1,nmax
756  cn=ci**(nn-n-1)
757  dnn=dfloat((2*n+1)*(2*nn+1))
758  dnn=dnn/dfloat( n*nn*(n+1)*(nn+1) )
759  rn=dsqrt(dnn)
760  cal(n,nn)=cn*rn
761  5 CONTINUE
762  dcth0=ctp
763  dcth=ctp1
764  ph=phip1-phip
765  vv=(0d0,0d0)
766  vh=(0d0,0d0)
767  hv=(0d0,0d0)
768  hh=(0d0,0d0)
769  DO 500 m=0,nmax
770  m1=m+1
771  nmin=max(m,1)
772  CALL vigampl (dcth, nmax, m, dv1, dv2)
773  CALL vigampl (dcth0, nmax, m, dv01, dv02)
774  fc=2d0*dcos(m*ph)
775  fs=2d0*dsin(m*ph)
776  DO 400 nn=nmin,nmax
777  dv1nn=m*dv01(nn)
778  dv2nn=dv02(nn)
779  DO 400 n=nmin,nmax
780  dv1n=m*dv1(n)
781  dv2n=dv2(n)
782 
783  ct11=dcmplx(tr11(m1,n,nn),ti11(m1,n,nn))
784  ct22=dcmplx(tr22(m1,n,nn),ti22(m1,n,nn))
785 
786  IF (m.EQ.0) THEN
787 
788  cn=cal(n,nn)*dv2n*dv2nn
789 
790  vv=vv+cn*ct22
791  hh=hh+cn*ct11
792 
793  ELSE
794 
795  ct12=dcmplx(tr12(m1,n,nn),ti12(m1,n,nn))
796  ct21=dcmplx(tr21(m1,n,nn),ti21(m1,n,nn))
797 
798  cn1=cal(n,nn)*fc
799  cn2=cal(n,nn)*fs
800 
801  d11=dv1n*dv1nn
802  d12=dv1n*dv2nn
803  d21=dv2n*dv1nn
804  d22=dv2n*dv2nn
805 
806  vv=vv+(ct11*d11+ct21*d21
807  & +ct12*d12+ct22*d22)*cn1
808 
809  vh=vh+(ct11*d12+ct21*d22
810  & +ct12*d11+ct22*d21)*cn2
811 
812  hv=hv-(ct11*d21+ct21*d11
813  & +ct12*d22+ct22*d12)*cn2
814 
815  hh=hh+(ct11*d22+ct21*d12
816  & +ct12*d21+ct22*d11)*cn1
817  ENDIF
818  400 CONTINUE
819  500 CONTINUE
820  dk=2d0*pin/dlam
821  vv=vv/dk
822  vh=vh/dk
823  hv=hv/dk
824  hh=hh/dk
825  cvv=vv*r(1,1)+vh*r(2,1)
826  cvh=vv*r(1,2)+vh*r(2,2)
827  chv=hv*r(1,1)+hh*r(2,1)
828  chh=hv*r(1,2)+hh*r(2,2)
829  vv=r1(1,1)*cvv+r1(1,2)*chv
830  vh=r1(1,1)*cvh+r1(1,2)*chh
831  hv=r1(2,1)*cvv+r1(2,2)*chv
832  hh=r1(2,1)*cvh+r1(2,2)*chh
833 
834 C WRITE (6,1005) TL,TL1,PL,PL1,ALPHA,BETA
835 C WRITE (6,1006)
836 C PRINT 1101, VV
837 C PRINT 1102, VH
838 C PRINT 1103, HV
839 C PRINT 1104, HH
840 C 1101 FORMAT ('S11=',D11.5,' + i*',D11.5)
841 C 1102 FORMAT ('S12=',D11.5,' + i*',D11.5)
842 C 1103 FORMAT ('S21=',D11.5,' + i*',D11.5)
843 C 1104 FORMAT ('S22=',D11.5,' + i*',D11.5)
844 C 1005 FORMAT ('thet0=',F6.2,' thet=',F6.2,' phi0=',F6.2,
845 C & ' phi=',F6.2,' alpha=',F6.2,' beta=',F6.2)
846 C 1006 FORMAT ('AMPLITUDE MATRIX')
847  RETURN
848  END
849 
850 C*****************************************************************
851 C
852 C Calculation of the functions
853 C DV1(N)=dvig(0,m,n,arccos x)/sin(arccos x)
854 C and
855 C DV2(N)=[d/d(arccos x)] dvig(0,m,n,arccos x)
856 C 1.LE.N.LE.NMAX
857 C 0.LE.X.LE.1
858 
859  SUBROUTINE vigampl (X, NMAX, M, DV1, DV2)
860  include 'tm_tmd.par.f'
861  IMPLICIT REAL*8 (a-h,o-z)
862  real*8 dv1(npn6), dv2(npn6)
863  DO 1 n=1,nmax
864  dv1(n)=0d0
865  dv2(n)=0d0
866  1 CONTINUE
867  dx=dabs(x)
868  IF (dabs(1d0-dx).LE.1d-10) GO TO 100
869  a=1d0
870  qs=dsqrt(1d0-x*x)
871  qs1=1d0/qs
872  dsi=qs1
873  IF (m.NE.0) GO TO 20
874  d1=1d0
875  d2=x
876  DO 5 n=1,nmax
877  qn=dfloat(n)
878  qn1=dfloat(n+1)
879  qn2=dfloat(2*n+1)
880  d3=(qn2*x*d2-qn*d1)/qn1
881  der=qs1*(qn1*qn/qn2)*(-d1+d3)
882  dv1(n)=d2*dsi
883  dv2(n)=der
884  d1=d2
885  d2=d3
886  5 CONTINUE
887  RETURN
888  20 qmm=dfloat(m*m)
889  DO 25 i=1,m
890  i2=i*2
891  a=a*dsqrt(dfloat(i2-1)/dfloat(i2))*qs
892  25 CONTINUE
893  d1=0d0
894  d2=a
895  DO 30 n=m,nmax
896  qn=dfloat(n)
897  qn2=dfloat(2*n+1)
898  qn1=dfloat(n+1)
899  qnm=dsqrt(qn*qn-qmm)
900  qnm1=dsqrt(qn1*qn1-qmm)
901  d3=(qn2*x*d2-qnm*d1)/qnm1
902  der=qs1*(-qn1*qnm*d1+qn*qnm1*d3)/qn2
903  dv1(n)=d2*dsi
904  dv2(n)=der
905  d1=d2
906  d2=d3
907  30 CONTINUE
908  RETURN
909  100 IF (m.NE.1) RETURN
910  DO 110 n=1,nmax
911  dn=dfloat(n*(n+1))
912  dn=0.5d0*dsqrt(dn)
913  IF (x.LT.0d0) dn=dn*(-1)**(n+1)
914  dv1(n)=dn
915  IF (x.LT.0d0) dn=-dn
916  dv2(n)=dn
917  110 CONTINUE
918  RETURN
919  END
920 
921 C**********************************************************************
922 
923  SUBROUTINE const (NGAUSS,NMAX,MMAX,P,X,W,AN,ANN,S,SS,NP,EPS)
924  IMPLICIT REAL*8 (a-h,o-z)
925  include 'tm_tmd.par.f'
926  real*8 x(npng2),w(npng2),x1(npng1),w1(npng1),
927  * x2(npng1),w2(npng1),
928  * s(npng2),ss(npng2),
929  * an(npn1),ann(npn1,npn1),dd(npn1)
930 
931  DO 10 n=1,nmax
932  nn=n*(n+1)
933  an(n)=dfloat(nn)
934  d=dsqrt(dfloat(2*n+1)/dfloat(nn))
935  dd(n)=d
936  DO 10 n1=1,n
937  ddd=d*dd(n1)*0.5d0
938  ann(n,n1)=ddd
939  ann(n1,n)=ddd
940  10 CONTINUE
941  ng=2*ngauss
942  IF (np.EQ.-2) GO TO 11
943  CALL gauss(ng,0,0,x,w)
944  GO TO 19
945  11 ng1=dfloat(ngauss)/2d0
946  ng2=ngauss-ng1
947  xx=-dcos(datan(eps))
948  CALL gauss(ng1,0,0,x1,w1)
949  CALL gauss(ng2,0,0,x2,w2)
950  DO 12 i=1,ng1
951  w(i)=0.5d0*(xx+1d0)*w1(i)
952  x(i)=0.5d0*(xx+1d0)*x1(i)+0.5d0*(xx-1d0)
953  12 CONTINUE
954  DO 14 i=1,ng2
955  w(i+ng1)=-0.5d0*xx*w2(i)
956  x(i+ng1)=-0.5d0*xx*x2(i)+0.5d0*xx
957  14 CONTINUE
958  DO 16 i=1,ngauss
959  w(ng-i+1)=w(i)
960  x(ng-i+1)=-x(i)
961  16 CONTINUE
962  19 DO 20 i=1,ngauss
963  y=x(i)
964  y=1d0/(1d0-y*y)
965  ss(i)=y
966  ss(ng-i+1)=y
967  y=dsqrt(y)
968  s(i)=y
969  s(ng-i+1)=y
970  20 CONTINUE
971  RETURN
972  END
973 
974 C**********************************************************************
975 
976  SUBROUTINE vary (LAM,MRR,MRI,A,EPS,NP,NGAUSS,X,P,PPI,PIR,PII,
977  * R,DR,DDR,DRR,DRI,NMAX)
978  include 'tm_tmd.par.f'
979  IMPLICIT REAL*8 (a-h,o-z)
980  real*8 x(npng2),r(npng2),dr(npng2),mrr,mri,lam,
981  * z(npng2),zr(npng2),zi(npng2),
982  * j(npng2,npn1),y(npng2,npn1),jr(npng2,npn1),
983  * ji(npng2,npn1),dj(npng2,npn1),
984  * djr(npng2,npn1),dji(npng2,npn1),ddr(npng2),
985  * drr(npng2),dri(npng2),
986  * dy(npng2,npn1)
987  COMMON /cbess/ j,y,jr,ji,dj,dy,djr,dji
988  ng=ngauss*2
989  IF (np.GT.0) CALL rsp2(x,ng,a,eps,np,r,dr)
990  IF (np.EQ.-1) CALL rsp1(x,ng,ngauss,a,eps,np,r,dr)
991  IF (np.EQ.-2) CALL rsp3(x,ng,ngauss,a,eps,r,dr)
992  IF (np.EQ.-3) CALL rsp4(x,ng,a,r,dr)
993  pi=p*2d0/lam
994  ppi=pi*pi
995  pir=ppi*mrr
996  pii=ppi*mri
997  v=1d0/(mrr*mrr+mri*mri)
998  prr=mrr*v
999  pri=-mri*v
1000  ta=0d0
1001  DO 10 i=1,ng
1002  vv=dsqrt(r(i))
1003  v=vv*pi
1004  ta=max(ta,v)
1005  vv=1d0/v
1006  ddr(i)=vv
1007  drr(i)=prr*vv
1008  dri(i)=pri*vv
1009  v1=v*mrr
1010  v2=v*mri
1011  z(i)=v
1012  zr(i)=v1
1013  zi(i)=v2
1014  10 CONTINUE
1015 C IF (NMAX.GT.NPN1) PRINT 9000,NMAX,NPN1
1016  IF (nmax.GT.npn1) stop
1017 C 9000 FORMAT(' NMAX = ',I2,', i.e., greater than ',I3)
1018  tb=ta*dsqrt(mrr*mrr+mri*mri)
1019  tb=dmax1(tb,dfloat(nmax))
1020  nnmax1=1.2d0*dsqrt(dmax1(ta,dfloat(nmax)))+3d0
1021  nnmax2=(tb+4d0*(tb**0.33333d0)+1.2d0*dsqrt(tb))
1022  nnmax2=nnmax2-nmax+5
1023  CALL bess(z,zr,zi,ng,nmax,nnmax1,nnmax2)
1024  RETURN
1025  END
1026 
1027 C**********************************************************************
1028 
1029  SUBROUTINE rsp1 (X,NG,NGAUSS,REV,EPS,NP,R,DR)
1030  IMPLICIT REAL*8 (a-h,o-z)
1031  real*8 x(ng),r(ng),dr(ng)
1032  a=rev*eps**(1d0/3d0)
1033  aa=a*a
1034  ee=eps*eps
1035  ee1=ee-1d0
1036  DO 50 i=1,ngauss
1037  c=x(i)
1038  cc=c*c
1039  ss=1d0-cc
1040  s=dsqrt(ss)
1041  rr=1d0/(ss+ee*cc)
1042  r(i)=aa*rr
1043  r(ng-i+1)=r(i)
1044  dr(i)=rr*c*s*ee1
1045  dr(ng-i+1)=-dr(i)
1046  50 CONTINUE
1047  RETURN
1048  END
1049 
1050 C**********************************************************************
1051 
1052  SUBROUTINE rsp2 (X,NG,REV,EPS,N,R,DR)
1053  IMPLICIT REAL*8 (a-h,o-z)
1054  real*8 x(ng),r(ng),dr(ng)
1055  dnp=dfloat(n)
1056  dn=dnp*dnp
1057  dn4=dn*4d0
1058  ep=eps*eps
1059  a=1d0+1.5d0*ep*(dn4-2d0)/(dn4-1d0)
1060  i=(dnp+0.1d0)*0.5d0
1061  i=2*i
1062  IF (i.EQ.n) a=a-3d0*eps*(1d0+0.25d0*ep)/
1063  * (dn-1d0)-0.25d0*ep*eps/(9d0*dn-1d0)
1064  r0=rev*a**(-1d0/3d0)
1065  DO 50 i=1,ng
1066  xi=dacos(x(i))*dnp
1067  ri=r0*(1d0+eps*dcos(xi))
1068  r(i)=ri*ri
1069  dr(i)=-r0*eps*dnp*dsin(xi)/ri
1070 c WRITE (6,*) I,R(I),DR(I)
1071  50 CONTINUE
1072  RETURN
1073  END
1074 
1075 C**********************************************************************
1076 
1077  SUBROUTINE rsp3 (X,NG,NGAUSS,REV,EPS,R,DR)
1078  IMPLICIT REAL*8 (a-h,o-z)
1079  real*8 x(ng),r(ng),dr(ng)
1080  h=rev*( (2d0/(3d0*eps*eps))**(1d0/3d0) )
1081  a=h*eps
1082  DO 50 i=1,ngauss
1083  co=-x(i)
1084  si=dsqrt(1d0-co*co)
1085  IF (si/co.GT.a/h) GO TO 20
1086  rad=h/co
1087  rthet=h*si/(co*co)
1088  GO TO 30
1089  20 rad=a/si
1090  rthet=-a*co/(si*si)
1091  30 r(i)=rad*rad
1092  r(ng-i+1)=r(i)
1093  dr(i)=-rthet/rad
1094  dr(ng-i+1)=-dr(i)
1095  50 CONTINUE
1096  RETURN
1097  END
1098 
1099 C**********************************************************************
1100 C *
1101 C Calculation of the functions R(I)=r(y)**2 and *
1102 C DR(I)=((d/dy)r(y))/r(y) for a distorted *
1103 C droplet specified by the parameters r_ev (equal-volume-sphere *
1104 C radius) and c_n (Chebyshev expansion coefficients) *
1105 C Y(I)=arccos(X(I)) *
1106 C 1.LE.I.LE.NG *
1107 C X - arguments *
1108 C *
1109 C**********************************************************************
1110 
1111  SUBROUTINE rsp4 (X,NG,REV,R,DR)
1112  parameter(nc=10)
1113  IMPLICIT REAL*8 (a-h,o-z)
1114  real*8 x(ng),r(ng),dr(ng),c(0:nc)
1115  COMMON /cdrop/ c,r0v
1116  r0=rev*r0v
1117  DO i=1,ng
1118  xi=dacos(x(i))
1119  ri=1d0+c(0)
1120  dri=0d0
1121  DO n=1,nc
1122  xin=xi*n
1123  ri=ri+c(n)*dcos(xin)
1124  dri=dri-c(n)*n*dsin(xin)
1125  ENDDO
1126  ri=ri*r0
1127  dri=dri*r0
1128  r(i)=ri*ri
1129  dr(i)=dri/ri
1130 c WRITE (6,*) I,R(I),DR(I)
1131  ENDDO
1132  RETURN
1133  END
1134 
1135 C*********************************************************************
1136 
1137  SUBROUTINE bess (X,XR,XI,NG,NMAX,NNMAX1,NNMAX2)
1138  include 'tm_tmd.par.f'
1139  IMPLICIT REAL*8 (a-h,o-z)
1140  real*8 x(ng),xr(ng),xi(ng),
1141  * j(npng2,npn1),y(npng2,npn1),jr(npng2,npn1),
1142  * ji(npng2,npn1),dj(npng2,npn1),dy(npng2,npn1),
1143  * djr(npng2,npn1),dji(npng2,npn1),
1144  * aj(npn1),ay(npn1),ajr(npn1),aji(npn1),
1145  * adj(npn1),ady(npn1),adjr(npn1),
1146  * adji(npn1)
1147  COMMON /cbess/ j,y,jr,ji,dj,dy,djr,dji
1148 
1149  DO 10 i=1,ng
1150  xx=x(i)
1151  CALL rjb(xx,aj,adj,nmax,nnmax1)
1152  CALL ryb(xx,ay,ady,nmax)
1153  yr=xr(i)
1154  yi=xi(i)
1155  CALL cjb(yr,yi,ajr,aji,adjr,adji,nmax,nnmax2)
1156  DO 10 n=1,nmax
1157  j(i,n)=aj(n)
1158  y(i,n)=ay(n)
1159  jr(i,n)=ajr(n)
1160  ji(i,n)=aji(n)
1161  dj(i,n)=adj(n)
1162  dy(i,n)=ady(n)
1163  djr(i,n)=adjr(n)
1164  dji(i,n)=adji(n)
1165  10 CONTINUE
1166  RETURN
1167  END
1168 
1169 C**********************************************************************
1170 
1171  SUBROUTINE rjb(X,Y,U,NMAX,NNMAX)
1172  IMPLICIT REAL*8 (a-h,o-z)
1173  real*8 y(nmax),u(nmax),z(800)
1174  l=nmax+nnmax
1175  xx=1d0/x
1176  z(l)=1d0/(dfloat(2*l+1)*xx)
1177  l1=l-1
1178  DO 5 i=1,l1
1179  i1=l-i
1180  z(i1)=1d0/(dfloat(2*i1+1)*xx-z(i1+1))
1181  5 CONTINUE
1182  z0=1d0/(xx-z(1))
1183  y0=z0*dcos(x)*xx
1184  y1=y0*z(1)
1185  u(1)=y0-y1*xx
1186  y(1)=y1
1187  DO 10 i=2,nmax
1188  yi1=y(i-1)
1189  yi=yi1*z(i)
1190  u(i)=yi1-dfloat(i)*yi*xx
1191  y(i)=yi
1192  10 CONTINUE
1193  RETURN
1194  END
1195 
1196 C**********************************************************************
1197 
1198  SUBROUTINE ryb(X,Y,V,NMAX)
1199  IMPLICIT REAL*8 (a-h,o-z)
1200  real*8 y(nmax),v(nmax)
1201  c=dcos(x)
1202  s=dsin(x)
1203  x1=1d0/x
1204  x2=x1*x1
1205  x3=x2*x1
1206  y1=-c*x2-s*x1
1207  y(1)=y1
1208  y(2)=(-3d0*x3+x1)*c-3d0*x2*s
1209  nmax1=nmax-1
1210  DO 5 i=2,nmax1
1211  5 y(i+1)=dfloat(2*i+1)*x1*y(i)-y(i-1)
1212  v(1)=-x1*(c+y1)
1213  DO 10 i=2,nmax
1214  10 v(i)=y(i-1)-dfloat(i)*x1*y(i)
1215  RETURN
1216  END
1217 
1218 C**********************************************************************
1219 C *
1220 C CALCULATION OF SPHERICAL BESSEL FUNCTIONS OF THE FIRST KIND *
1221 C J=JR+I*JI OF COMPLEX ARGUMENT X=XR+I*XI OF ORDERS FROM 1 TO NMAX *
1222 C BY USING BACKWARD RECURSION. PARAMETR NNMAX DETERMINES NUMERICAL *
1223 C ACCURACY. U=UR+I*UI - FUNCTION (1/X)(D/DX)(X*J(X)) *
1224 C *
1225 C**********************************************************************
1226 
1227  SUBROUTINE cjb (XR,XI,YR,YI,UR,UI,NMAX,NNMAX)
1228  include 'tm_tmd.par.f'
1229  IMPLICIT REAL*8 (a-h,o-z)
1230  real*8 yr(nmax),yi(nmax),ur(nmax),ui(nmax)
1231  real*8 cyr(npn1),cyi(npn1),czr(1200),czi(1200),
1232  * cur(npn1),cui(npn1)
1233  l=nmax+nnmax
1234  xrxi=1d0/(xr*xr+xi*xi)
1235  cxxr=xr*xrxi
1236  cxxi=-xi*xrxi
1237  qf=1d0/dfloat(2*l+1)
1238  czr(l)=xr*qf
1239  czi(l)=xi*qf
1240  l1=l-1
1241  DO i=1,l1
1242  i1=l-i
1243  qf=dfloat(2*i1+1)
1244  ar=qf*cxxr-czr(i1+1)
1245  ai=qf*cxxi-czi(i1+1)
1246  ari=1d0/(ar*ar+ai*ai)
1247  czr(i1)=ar*ari
1248  czi(i1)=-ai*ari
1249  ENDDO
1250  ar=cxxr-czr(1)
1251  ai=cxxi-czi(1)
1252  ari=1d0/(ar*ar+ai*ai)
1253  cz0r=ar*ari
1254  cz0i=-ai*ari
1255  cr=dcos(xr)*dcosh(xi)
1256  ci=-dsin(xr)*dsinh(xi)
1257  ar=cz0r*cr-cz0i*ci
1258  ai=cz0i*cr+cz0r*ci
1259  cy0r=ar*cxxr-ai*cxxi
1260  cy0i=ai*cxxr+ar*cxxi
1261  cy1r=cy0r*czr(1)-cy0i*czi(1)
1262  cy1i=cy0i*czr(1)+cy0r*czi(1)
1263  ar=cy1r*cxxr-cy1i*cxxi
1264  ai=cy1i*cxxr+cy1r*cxxi
1265  cu1r=cy0r-ar
1266  cu1i=cy0i-ai
1267  cyr(1)=cy1r
1268  cyi(1)=cy1i
1269  cur(1)=cu1r
1270  cui(1)=cu1i
1271  yr(1)=cy1r
1272  yi(1)=cy1i
1273  ur(1)=cu1r
1274  ui(1)=cu1i
1275  DO i=2,nmax
1276  qi=dfloat(i)
1277  cyi1r=cyr(i-1)
1278  cyi1i=cyi(i-1)
1279  cyir=cyi1r*czr(i)-cyi1i*czi(i)
1280  cyii=cyi1i*czr(i)+cyi1r*czi(i)
1281  ar=cyir*cxxr-cyii*cxxi
1282  ai=cyii*cxxr+cyir*cxxi
1283  cuir=cyi1r-qi*ar
1284  cuii=cyi1i-qi*ai
1285  cyr(i)=cyir
1286  cyi(i)=cyii
1287  cur(i)=cuir
1288  cui(i)=cuii
1289  yr(i)=cyir
1290  yi(i)=cyii
1291  ur(i)=cuir
1292  ui(i)=cuii
1293  ENDDO
1294  RETURN
1295  END
1296 
1297 C**********************************************************************
1298 
1299  SUBROUTINE tmatr0 (NGAUSS,X,W,AN,ANN,S,SS,PPI,PIR,PII,R,DR,DDR,
1300  * DRR,DRI,NMAX,NCHECK)
1301  include 'tm_tmd.par.f'
1302  IMPLICIT REAL*8 (a-h,o-z)
1303  real*8 x(npng2),w(npng2),an(npn1),s(npng2),ss(npng2),
1304  * r(npng2),dr(npng2),sig(npn2),
1305  * j(npng2,npn1),y(npng2,npn1),
1306  * jr(npng2,npn1),ji(npng2,npn1),dj(npng2,npn1),
1307  * dy(npng2,npn1),djr(npng2,npn1),
1308  * dji(npng2,npn1),ddr(npng2),drr(npng2),
1309  * d1(npng2,npn1),d2(npng2,npn1),
1310  * dri(npng2),ds(npng2),dss(npng2),rr(npng2),
1311  * dv1(npn1),dv2(npn1)
1312 
1313  real*8 r11(npn1,npn1),r12(npn1,npn1),
1314  * r21(npn1,npn1),r22(npn1,npn1),
1315  * i11(npn1,npn1),i12(npn1,npn1),
1316  * i21(npn1,npn1),i22(npn1,npn1),
1317  * rg11(npn1,npn1),rg12(npn1,npn1),
1318  * rg21(npn1,npn1),rg22(npn1,npn1),
1319  * ig11(npn1,npn1),ig12(npn1,npn1),
1320  * ig21(npn1,npn1),ig22(npn1,npn1),
1321  * ann(npn1,npn1),
1322  * qr(npn2,npn2),qi(npn2,npn2),
1323  * rgqr(npn2,npn2),rgqi(npn2,npn2),
1324  * tqr(npn2,npn2),tqi(npn2,npn2),
1325  * trgqr(npn2,npn2),trgqi(npn2,npn2)
1326  real*8 tr1(npn2,npn2),ti1(npn2,npn2)
1327  COMMON /tmat99/
1328  & r11,r12,r21,r22,i11,i12,i21,i22,rg11,rg12,rg21,rg22,
1329  & ig11,ig12,ig21,ig22
1330  COMMON /cbess/ j,y,jr,ji,dj,dy,djr,dji
1331  COMMON /ct/ tr1,ti1
1332  COMMON /ctt/ qr,qi,rgqr,rgqi
1333  mm1=1
1334  nnmax=nmax+nmax
1335  ng=2*ngauss
1336  ngss=ng
1337  factor=1d0
1338  IF (ncheck.EQ.1) THEN
1339  ngss=ngauss
1340  factor=2d0
1341  ELSE
1342  CONTINUE
1343  ENDIF
1344  si=1d0
1345  DO 5 n=1,nnmax
1346  si=-si
1347  sig(n)=si
1348  5 CONTINUE
1349  20 DO 25 i=1,ngauss
1350  i1=ngauss+i
1351  i2=ngauss-i+1
1352  CALL vig ( x(i1), nmax, 0, dv1, dv2)
1353  DO 25 n=1,nmax
1354  si=sig(n)
1355  dd1=dv1(n)
1356  dd2=dv2(n)
1357  d1(i1,n)=dd1
1358  d2(i1,n)=dd2
1359  d1(i2,n)=dd1*si
1360  d2(i2,n)=-dd2*si
1361  25 CONTINUE
1362  30 DO 40 i=1,ngss
1363  rr(i)=w(i)*r(i)
1364  40 CONTINUE
1365 
1366  DO 300 n1=mm1,nmax
1367  an1=an(n1)
1368  DO 300 n2=mm1,nmax
1369  an2=an(n2)
1370  ar12=0d0
1371  ar21=0d0
1372  ai12=0d0
1373  ai21=0d0
1374  gr12=0d0
1375  gr21=0d0
1376  gi12=0d0
1377  gi21=0d0
1378  IF (ncheck.EQ.1.AND.sig(n1+n2).LT.0d0) GO TO 205
1379 
1380  DO 200 i=1,ngss
1381  d1n1=d1(i,n1)
1382  d2n1=d2(i,n1)
1383  d1n2=d1(i,n2)
1384  d2n2=d2(i,n2)
1385  a12=d1n1*d2n2
1386  a21=d2n1*d1n2
1387  a22=d2n1*d2n2
1388  aa1=a12+a21
1389 
1390  qj1=j(i,n1)
1391  qy1=y(i,n1)
1392  qjr2=jr(i,n2)
1393  qji2=ji(i,n2)
1394  qdjr2=djr(i,n2)
1395  qdji2=dji(i,n2)
1396  qdj1=dj(i,n1)
1397  qdy1=dy(i,n1)
1398 
1399  c1r=qjr2*qj1
1400  c1i=qji2*qj1
1401  b1r=c1r-qji2*qy1
1402  b1i=c1i+qjr2*qy1
1403 
1404  c2r=qjr2*qdj1
1405  c2i=qji2*qdj1
1406  b2r=c2r-qji2*qdy1
1407  b2i=c2i+qjr2*qdy1
1408 
1409  ddri=ddr(i)
1410  c3r=ddri*c1r
1411  c3i=ddri*c1i
1412  b3r=ddri*b1r
1413  b3i=ddri*b1i
1414 
1415  c4r=qdjr2*qj1
1416  c4i=qdji2*qj1
1417  b4r=c4r-qdji2*qy1
1418  b4i=c4i+qdjr2*qy1
1419 
1420  drri=drr(i)
1421  drii=dri(i)
1422  c5r=c1r*drri-c1i*drii
1423  c5i=c1i*drri+c1r*drii
1424  b5r=b1r*drri-b1i*drii
1425  b5i=b1i*drri+b1r*drii
1426 
1427  uri=dr(i)
1428  rri=rr(i)
1429 
1430  f1=rri*a22
1431  f2=rri*uri*an1*a12
1432  ar12=ar12+f1*b2r+f2*b3r
1433  ai12=ai12+f1*b2i+f2*b3i
1434  gr12=gr12+f1*c2r+f2*c3r
1435  gi12=gi12+f1*c2i+f2*c3i
1436 
1437  f2=rri*uri*an2*a21
1438  ar21=ar21+f1*b4r+f2*b5r
1439  ai21=ai21+f1*b4i+f2*b5i
1440  gr21=gr21+f1*c4r+f2*c5r
1441  gi21=gi21+f1*c4i+f2*c5i
1442  200 CONTINUE
1443 
1444  205 an12=ann(n1,n2)*factor
1445  r12(n1,n2)=ar12*an12
1446  r21(n1,n2)=ar21*an12
1447  i12(n1,n2)=ai12*an12
1448  i21(n1,n2)=ai21*an12
1449  rg12(n1,n2)=gr12*an12
1450  rg21(n1,n2)=gr21*an12
1451  ig12(n1,n2)=gi12*an12
1452  ig21(n1,n2)=gi21*an12
1453  300 CONTINUE
1454 
1455  tpir=pir
1456  tpii=pii
1457  tppi=ppi
1458 
1459  nm=nmax
1460  DO 310 n1=mm1,nmax
1461  k1=n1-mm1+1
1462  kk1=k1+nm
1463  DO 310 n2=mm1,nmax
1464  k2=n2-mm1+1
1465  kk2=k2+nm
1466 
1467  tar12= i12(n1,n2)
1468  tai12=-r12(n1,n2)
1469  tgr12= ig12(n1,n2)
1470  tgi12=-rg12(n1,n2)
1471 
1472  tar21=-i21(n1,n2)
1473  tai21= r21(n1,n2)
1474  tgr21=-ig21(n1,n2)
1475  tgi21= rg21(n1,n2)
1476 
1477  tqr(k1,k2)=tpir*tar21-tpii*tai21+tppi*tar12
1478  tqi(k1,k2)=tpir*tai21+tpii*tar21+tppi*tai12
1479  trgqr(k1,k2)=tpir*tgr21-tpii*tgi21+tppi*tgr12
1480  trgqi(k1,k2)=tpir*tgi21+tpii*tgr21+tppi*tgi12
1481 
1482  tqr(k1,kk2)=0d0
1483  tqi(k1,kk2)=0d0
1484  trgqr(k1,kk2)=0d0
1485  trgqi(k1,kk2)=0d0
1486 
1487  tqr(kk1,k2)=0d0
1488  tqi(kk1,k2)=0d0
1489  trgqr(kk1,k2)=0d0
1490  trgqi(kk1,k2)=0d0
1491 
1492  tqr(kk1,kk2)=tpir*tar12-tpii*tai12+tppi*tar21
1493  tqi(kk1,kk2)=tpir*tai12+tpii*tar12+tppi*tai21
1494  trgqr(kk1,kk2)=tpir*tgr12-tpii*tgi12+tppi*tgr21
1495  trgqi(kk1,kk2)=tpir*tgi12+tpii*tgr12+tppi*tgi21
1496  310 CONTINUE
1497 
1498  nnmax=2*nm
1499  DO 320 n1=1,nnmax
1500  DO 320 n2=1,nnmax
1501  qr(n1,n2)=tqr(n1,n2)
1502  qi(n1,n2)=tqi(n1,n2)
1503  rgqr(n1,n2)=trgqr(n1,n2)
1504  rgqi(n1,n2)=trgqi(n1,n2)
1505  320 CONTINUE
1506  CALL tt(nmax,ncheck)
1507  RETURN
1508  END
1509 
1510 C**********************************************************************
1511 
1512  SUBROUTINE tmatr (M,NGAUSS,X,W,AN,ANN,S,SS,PPI,PIR,PII,R,DR,DDR,
1513  * DRR,DRI,NMAX,NCHECK)
1514  include 'tm_tmd.par.f'
1515  IMPLICIT REAL*8 (a-h,o-z)
1516  real*8 x(npng2),w(npng2),an(npn1),s(npng2),ss(npng2),
1517  * r(npng2),dr(npng2),sig(npn2),
1518  * j(npng2,npn1),y(npng2,npn1),
1519  * jr(npng2,npn1),ji(npng2,npn1),dj(npng2,npn1),
1520  * dy(npng2,npn1),djr(npng2,npn1),
1521  * dji(npng2,npn1),ddr(npng2),drr(npng2),
1522  * d1(npng2,npn1),d2(npng2,npn1),
1523  * dri(npng2),ds(npng2),dss(npng2),rr(npng2),
1524  * dv1(npn1),dv2(npn1)
1525 
1526  real*8 r11(npn1,npn1),r12(npn1,npn1),
1527  * r21(npn1,npn1),r22(npn1,npn1),
1528  * i11(npn1,npn1),i12(npn1,npn1),
1529  * i21(npn1,npn1),i22(npn1,npn1),
1530  * rg11(npn1,npn1),rg12(npn1,npn1),
1531  * rg21(npn1,npn1),rg22(npn1,npn1),
1532  * ig11(npn1,npn1),ig12(npn1,npn1),
1533  * ig21(npn1,npn1),ig22(npn1,npn1),
1534  * ann(npn1,npn1),
1535  * qr(npn2,npn2),qi(npn2,npn2),
1536  * rgqr(npn2,npn2),rgqi(npn2,npn2),
1537  * tqr(npn2,npn2),tqi(npn2,npn2),
1538  * trgqr(npn2,npn2),trgqi(npn2,npn2)
1539  real*8 tr1(npn2,npn2),ti1(npn2,npn2)
1540  COMMON /tmat99/
1541  & r11,r12,r21,r22,i11,i12,i21,i22,rg11,rg12,rg21,rg22,
1542  & ig11,ig12,ig21,ig22
1543  COMMON /cbess/ j,y,jr,ji,dj,dy,djr,dji
1544  COMMON /ct/ tr1,ti1
1545  COMMON /ctt/ qr,qi,rgqr,rgqi
1546  mm1=m
1547  qm=dfloat(m)
1548  qmm=qm*qm
1549  ng=2*ngauss
1550  ngss=ng
1551  factor=1d0
1552  IF (ncheck.EQ.1) THEN
1553  ngss=ngauss
1554  factor=2d0
1555  ELSE
1556  CONTINUE
1557  ENDIF
1558  si=1d0
1559  nm=nmax+nmax
1560  DO 5 n=1,nm
1561  si=-si
1562  sig(n)=si
1563  5 CONTINUE
1564  20 DO 25 i=1,ngauss
1565  i1=ngauss+i
1566  i2=ngauss-i+1
1567  CALL vig (x(i1),nmax,m,dv1,dv2)
1568  DO 25 n=1,nmax
1569  si=sig(n)
1570  dd1=dv1(n)
1571  dd2=dv2(n)
1572  d1(i1,n)=dd1
1573  d2(i1,n)=dd2
1574  d1(i2,n)=dd1*si
1575  d2(i2,n)=-dd2*si
1576  25 CONTINUE
1577  30 DO 40 i=1,ngss
1578  wr=w(i)*r(i)
1579  ds(i)=s(i)*qm*wr
1580  dss(i)=ss(i)*qmm
1581  rr(i)=wr
1582  40 CONTINUE
1583 
1584  DO 300 n1=mm1,nmax
1585  an1=an(n1)
1586  DO 300 n2=mm1,nmax
1587  an2=an(n2)
1588  ar11=0d0
1589  ar12=0d0
1590  ar21=0d0
1591  ar22=0d0
1592  ai11=0d0
1593  ai12=0d0
1594  ai21=0d0
1595  ai22=0d0
1596  gr11=0d0
1597  gr12=0d0
1598  gr21=0d0
1599  gr22=0d0
1600  gi11=0d0
1601  gi12=0d0
1602  gi21=0d0
1603  gi22=0d0
1604  si=sig(n1+n2)
1605 
1606  DO 200 i=1,ngss
1607  d1n1=d1(i,n1)
1608  d2n1=d2(i,n1)
1609  d1n2=d1(i,n2)
1610  d2n2=d2(i,n2)
1611  a11=d1n1*d1n2
1612  a12=d1n1*d2n2
1613  a21=d2n1*d1n2
1614  a22=d2n1*d2n2
1615  aa1=a12+a21
1616  aa2=a11*dss(i)+a22
1617  qj1=j(i,n1)
1618  qy1=y(i,n1)
1619  qjr2=jr(i,n2)
1620  qji2=ji(i,n2)
1621  qdjr2=djr(i,n2)
1622  qdji2=dji(i,n2)
1623  qdj1=dj(i,n1)
1624  qdy1=dy(i,n1)
1625 
1626  c1r=qjr2*qj1
1627  c1i=qji2*qj1
1628  b1r=c1r-qji2*qy1
1629  b1i=c1i+qjr2*qy1
1630 
1631  c2r=qjr2*qdj1
1632  c2i=qji2*qdj1
1633  b2r=c2r-qji2*qdy1
1634  b2i=c2i+qjr2*qdy1
1635 
1636  ddri=ddr(i)
1637  c3r=ddri*c1r
1638  c3i=ddri*c1i
1639  b3r=ddri*b1r
1640  b3i=ddri*b1i
1641 
1642  c4r=qdjr2*qj1
1643  c4i=qdji2*qj1
1644  b4r=c4r-qdji2*qy1
1645  b4i=c4i+qdjr2*qy1
1646 
1647  drri=drr(i)
1648  drii=dri(i)
1649  c5r=c1r*drri-c1i*drii
1650  c5i=c1i*drri+c1r*drii
1651  b5r=b1r*drri-b1i*drii
1652  b5i=b1i*drri+b1r*drii
1653 
1654  c6r=qdjr2*qdj1
1655  c6i=qdji2*qdj1
1656  b6r=c6r-qdji2*qdy1
1657  b6i=c6i+qdjr2*qdy1
1658 
1659  c7r=c4r*ddri
1660  c7i=c4i*ddri
1661  b7r=b4r*ddri
1662  b7i=b4i*ddri
1663 
1664  c8r=c2r*drri-c2i*drii
1665  c8i=c2i*drri+c2r*drii
1666  b8r=b2r*drri-b2i*drii
1667  b8i=b2i*drri+b2r*drii
1668 
1669  uri=dr(i)
1670  dsi=ds(i)
1671  dssi=dss(i)
1672  rri=rr(i)
1673 
1674  IF (ncheck.EQ.1.AND.si.GT.0d0) GO TO 150
1675 
1676  e1=dsi*aa1
1677  ar11=ar11+e1*b1r
1678  ai11=ai11+e1*b1i
1679  gr11=gr11+e1*c1r
1680  gi11=gi11+e1*c1i
1681  IF (ncheck.EQ.1) GO TO 160
1682 
1683  150 f1=rri*aa2
1684  f2=rri*uri*an1*a12
1685  ar12=ar12+f1*b2r+f2*b3r
1686  ai12=ai12+f1*b2i+f2*b3i
1687  gr12=gr12+f1*c2r+f2*c3r
1688  gi12=gi12+f1*c2i+f2*c3i
1689 
1690  f2=rri*uri*an2*a21
1691  ar21=ar21+f1*b4r+f2*b5r
1692  ai21=ai21+f1*b4i+f2*b5i
1693  gr21=gr21+f1*c4r+f2*c5r
1694  gi21=gi21+f1*c4i+f2*c5i
1695  IF (ncheck.EQ.1) GO TO 200
1696 
1697  160 e2=dsi*uri*a11
1698  e3=e2*an2
1699  e2=e2*an1
1700  ar22=ar22+e1*b6r+e2*b7r+e3*b8r
1701  ai22=ai22+e1*b6i+e2*b7i+e3*b8i
1702  gr22=gr22+e1*c6r+e2*c7r+e3*c8r
1703  gi22=gi22+e1*c6i+e2*c7i+e3*c8i
1704  200 CONTINUE
1705  an12=ann(n1,n2)*factor
1706  r11(n1,n2)=ar11*an12
1707  r12(n1,n2)=ar12*an12
1708  r21(n1,n2)=ar21*an12
1709  r22(n1,n2)=ar22*an12
1710  i11(n1,n2)=ai11*an12
1711  i12(n1,n2)=ai12*an12
1712  i21(n1,n2)=ai21*an12
1713  i22(n1,n2)=ai22*an12
1714  rg11(n1,n2)=gr11*an12
1715  rg12(n1,n2)=gr12*an12
1716  rg21(n1,n2)=gr21*an12
1717  rg22(n1,n2)=gr22*an12
1718  ig11(n1,n2)=gi11*an12
1719  ig12(n1,n2)=gi12*an12
1720  ig21(n1,n2)=gi21*an12
1721  ig22(n1,n2)=gi22*an12
1722 
1723  300 CONTINUE
1724  tpir=pir
1725  tpii=pii
1726  tppi=ppi
1727  nm=nmax-mm1+1
1728  DO 310 n1=mm1,nmax
1729  k1=n1-mm1+1
1730  kk1=k1+nm
1731  DO 310 n2=mm1,nmax
1732  k2=n2-mm1+1
1733  kk2=k2+nm
1734 
1735  tar11=-r11(n1,n2)
1736  tai11=-i11(n1,n2)
1737  tgr11=-rg11(n1,n2)
1738  tgi11=-ig11(n1,n2)
1739 
1740  tar12= i12(n1,n2)
1741  tai12=-r12(n1,n2)
1742  tgr12= ig12(n1,n2)
1743  tgi12=-rg12(n1,n2)
1744 
1745  tar21=-i21(n1,n2)
1746  tai21= r21(n1,n2)
1747  tgr21=-ig21(n1,n2)
1748  tgi21= rg21(n1,n2)
1749 
1750  tar22=-r22(n1,n2)
1751  tai22=-i22(n1,n2)
1752  tgr22=-rg22(n1,n2)
1753  tgi22=-ig22(n1,n2)
1754 
1755  tqr(k1,k2)=tpir*tar21-tpii*tai21+tppi*tar12
1756  tqi(k1,k2)=tpir*tai21+tpii*tar21+tppi*tai12
1757  trgqr(k1,k2)=tpir*tgr21-tpii*tgi21+tppi*tgr12
1758  trgqi(k1,k2)=tpir*tgi21+tpii*tgr21+tppi*tgi12
1759 
1760  tqr(k1,kk2)=tpir*tar11-tpii*tai11+tppi*tar22
1761  tqi(k1,kk2)=tpir*tai11+tpii*tar11+tppi*tai22
1762  trgqr(k1,kk2)=tpir*tgr11-tpii*tgi11+tppi*tgr22
1763  trgqi(k1,kk2)=tpir*tgi11+tpii*tgr11+tppi*tgi22
1764 
1765  tqr(kk1,k2)=tpir*tar22-tpii*tai22+tppi*tar11
1766  tqi(kk1,k2)=tpir*tai22+tpii*tar22+tppi*tai11
1767  trgqr(kk1,k2)=tpir*tgr22-tpii*tgi22+tppi*tgr11
1768  trgqi(kk1,k2)=tpir*tgi22+tpii*tgr22+tppi*tgi11
1769 
1770  tqr(kk1,kk2)=tpir*tar12-tpii*tai12+tppi*tar21
1771  tqi(kk1,kk2)=tpir*tai12+tpii*tar12+tppi*tai21
1772  trgqr(kk1,kk2)=tpir*tgr12-tpii*tgi12+tppi*tgr21
1773  trgqi(kk1,kk2)=tpir*tgi12+tpii*tgr12+tppi*tgi21
1774  310 CONTINUE
1775 
1776  nnmax=2*nm
1777  DO 320 n1=1,nnmax
1778  DO 320 n2=1,nnmax
1779  qr(n1,n2)=tqr(n1,n2)
1780  qi(n1,n2)=tqi(n1,n2)
1781  rgqr(n1,n2)=trgqr(n1,n2)
1782  rgqi(n1,n2)=trgqi(n1,n2)
1783  320 CONTINUE
1784 
1785  CALL tt(nm,ncheck)
1786 
1787  RETURN
1788  END
1789 
1790 C*****************************************************************
1791 
1792  SUBROUTINE vig (X, NMAX, M, DV1, DV2)
1793  include 'tm_tmd.par.f'
1794  IMPLICIT REAL*8 (a-h,o-z)
1795  real*8 dv1(npn1),dv2(npn1)
1796 
1797  a=1d0
1798  qs=dsqrt(1d0-x*x)
1799  qs1=1d0/qs
1800  DO n=1,nmax
1801  dv1(n)=0d0
1802  dv2(n)=0d0
1803  ENDDO
1804  IF (m.NE.0) GO TO 20
1805  d1=1d0
1806  d2=x
1807  DO n=1,nmax
1808  qn=dfloat(n)
1809  qn1=dfloat(n+1)
1810  qn2=dfloat(2*n+1)
1811  d3=(qn2*x*d2-qn*d1)/qn1
1812  der=qs1*(qn1*qn/qn2)*(-d1+d3)
1813  dv1(n)=d2
1814  dv2(n)=der
1815  d1=d2
1816  d2=d3
1817  ENDDO
1818  RETURN
1819  20 qmm=dfloat(m*m)
1820  DO i=1,m
1821  i2=i*2
1822  a=a*dsqrt(dfloat(i2-1)/dfloat(i2))*qs
1823  ENDDO
1824  d1=0d0
1825  d2=a
1826  DO n=m,nmax
1827  qn=dfloat(n)
1828  qn2=dfloat(2*n+1)
1829  qn1=dfloat(n+1)
1830  qnm=dsqrt(qn*qn-qmm)
1831  qnm1=dsqrt(qn1*qn1-qmm)
1832  d3=(qn2*x*d2-qnm*d1)/qnm1
1833  der=qs1*(-qn1*qnm*d1+qn*qnm1*d3)/qn2
1834  dv1(n)=d2
1835  dv2(n)=der
1836  d1=d2
1837  d2=d3
1838  ENDDO
1839  RETURN
1840  END
1841 
1842 C**********************************************************************
1843 C *
1844 C CALCULATION OF THE MATRIX T = - RG(Q) * (Q**(-1)) *
1845 C *
1846 C INPUT INFORTMATION IS IN COMMON /CTT/ *
1847 C OUTPUT INFORMATION IS IN COMMON /CT/ *
1848 C *
1849 C**********************************************************************
1850 
1851  SUBROUTINE tt(NMAX,NCHECK)
1852  include 'tm_tmd.par.f'
1853  IMPLICIT REAL*8 (a-h,o-z)
1854  real*8 f(npn2,npn2),b(npn2),work(npn2),
1855  * qr(npn2,npn2),qi(npn2,npn2),
1856  * rgqr(npn2,npn2),rgqi(npn2,npn2),
1857  * a(npn2,npn2),c(npn2,npn2),d(npn2,npn2),e(npn2,npn2)
1858  real*8 tr1(npn2,npn2),ti1(npn2,npn2)
1859  COMPLEX*16 ZQ(NPN2,NPN2),ZW(NPN2)
1860  INTEGER IPIV(NPN2),IPVT(NPN2)
1861  COMMON /CT/ TR1,TI1
1862  COMMON /ctt/ qr,qi,rgqr,rgqi
1863  ndim=npn2
1864  nnmax=2*nmax
1865 
1866 C Matrix inversion from LAPACK
1867 
1868  DO i=1,nnmax
1869  DO j=1,nnmax
1870  zq(i,j)=dcmplx(qr(i,j),qi(i,j))
1871  ENDDO
1872  ENDDO
1873  info=0
1874  CALL zgetrf(nnmax,nnmax,zq,npn2,ipiv,info)
1875 C IF (INFO.NE.0) WRITE (6,1100) INFO
1876  CALL zgetri(nnmax,zq,npn2,ipiv,zw,npn2,info)
1877 C IF (INFO.NE.0) WRITE (6,1100) INFO
1878 
1879 C 1100 FORMAT ('WARNING: info=', i2)
1880  DO i=1,nnmax
1881  DO j=1,nnmax
1882  tr=0d0
1883  ti=0d0
1884  DO k=1,nnmax
1885  arr=rgqr(i,k)
1886  ari=rgqi(i,k)
1887  ar=zq(k,j)
1888  ai=dimag(zq(k,j))
1889  tr=tr-arr*ar+ari*ai
1890  ti=ti-arr*ai-ari*ar
1891  ENDDO
1892  tr1(i,j)=tr
1893  ti1(i,j)=ti
1894  ENDDO
1895  ENDDO
1896  RETURN
1897  END
1898 
1899 C*****************************************************************
1900 
1901  SUBROUTINE sarea (D,RAT)
1902  IMPLICIT REAL*8 (a-h,o-z)
1903  IF (d.GE.1) GO TO 10
1904  e=dsqrt(1d0-d*d)
1905  r=0.5d0*(d**(2d0/3d0) + d**(-1d0/3d0)*dasin(e)/e)
1906  r=dsqrt(r)
1907  rat=1d0/r
1908  RETURN
1909  10 e=dsqrt(1d0-1d0/(d*d))
1910  r=0.25d0*(2d0*d**(2d0/3d0) + d**(-4d0/3d0)*dlog((1d0+e)/(1d0-e))
1911  & /e)
1912  r=dsqrt(r)
1913  rat=1d0/r
1914  return
1915  END
1916 
1917 c****************************************************************
1918 
1919  SUBROUTINE surfch (N,E,RAT)
1920  IMPLICIT REAL*8 (a-h,o-z)
1921  real*8 x(60),w(60)
1922  dn=dfloat(n)
1923  e2=e*e
1924  en=e*dn
1925  ng=60
1926  CALL gauss (ng,0,0,x,w)
1927  s=0d0
1928  v=0d0
1929  DO 10 i=1,ng
1930  xi=x(i)
1931  dx=dacos(xi)
1932  dxn=dn*dx
1933  ds=dsin(dx)
1934  dsn=dsin(dxn)
1935  dcn=dcos(dxn)
1936  a=1d0+e*dcn
1937  a2=a*a
1938  ens=en*dsn
1939  s=s+w(i)*a*dsqrt(a2+ens*ens)
1940  v=v+w(i)*(ds*a+xi*ens)*ds*a2
1941  10 CONTINUE
1942  rs=dsqrt(s*0.5d0)
1943  rv=(v*3d0/4d0)**(1d0/3d0)
1944  rat=rv/rs
1945  RETURN
1946  END
1947 
1948 C********************************************************************
1949 
1950  SUBROUTINE sareac (EPS,RAT)
1951  IMPLICIT REAL*8 (a-h,o-z)
1952  rat=(1.5d0/eps)**(1d0/3d0)
1953  rat=rat/dsqrt( (eps+2d0)/(2d0*eps) )
1954  RETURN
1955  END
1956 
1957 C**********************************************************************
1958 
1959  SUBROUTINE drop (RAT)
1960  parameter(nc=10, ng=60)
1961  IMPLICIT REAL*8 (a-h,o-z)
1962  real*8 x(ng),w(ng),c(0:nc)
1963  COMMON /cdrop/ c,r0v
1964  c(0)=-0.0481 d0
1965  c(1)= 0.0359 d0
1966  c(2)=-0.1263 d0
1967  c(3)= 0.0244 d0
1968  c(4)= 0.0091 d0
1969  c(5)=-0.0099 d0
1970  c(6)= 0.0015 d0
1971  c(7)= 0.0025 d0
1972  c(8)=-0.0016 d0
1973  c(9)=-0.0002 d0
1974  c(10)= 0.0010 d0
1975  CALL gauss (ng,0,0,x,w)
1976  s=0d0
1977  v=0d0
1978  DO i=1,ng
1979  xi=dacos(x(i))
1980  wi=w(i)
1981  ri=1d0+c(0)
1982  dri=0d0
1983  DO n=1,nc
1984  xin=xi*n
1985  ri=ri+c(n)*dcos(xin)
1986  dri=dri-c(n)*n*dsin(xin)
1987  ENDDO
1988  si=dsin(xi)
1989  ci=x(i)
1990  risi=ri*si
1991  s=s+wi*ri*dsqrt(ri*ri+dri*dri)
1992  v=v+wi*ri*risi*(risi-dri*ci)
1993  ENDDO
1994  rs=dsqrt(s*0.5d0)
1995  rv=(v*3d0*0.25d0)**(1d0/3d0)
1996  IF (dabs(rat-1d0).GT.1d-8) rat=rv/rs
1997  r0v=1d0/rv
1998 C WRITE (6,1000) R0V
1999  DO n=0,nc
2000 C WRITE (6,1001) N,C(N)
2001  ENDDO
2002 C 1000 FORMAT ('r_0/r_ev=',F7.4)
2003 C 1001 FORMAT ('c_',I2,'=',F7.4)
2004  RETURN
2005  END
2006 
2007 C**********************************************************************
2008 C CALCULATION OF POINTS AND WEIGHTS OF GAUSSIAN QUADRATURE *
2009 C FORMULA. IF IND1 = 0 - ON INTERVAL (-1,1), IF IND1 = 1 - ON *
2010 C INTERVAL (0,1). IF IND2 = 1 RESULTS ARE PRINTED. *
2011 C N - NUMBER OF POINTS *
2012 C Z - DIVISION POINTS *
2013 C W - WEIGHTS *
2014 **********************************************************************
2015 
2016  SUBROUTINE gauss (N,IND1,IND2,Z,W)
2017  IMPLICIT REAL*8 (a-h,p-z)
2018  real*8 z(n),w(n)
2019  a=1d0
2020  b=2d0
2021  c=3d0
2022  ind=mod(n,2)
2023  k=n/2+ind
2024  f=dfloat(n)
2025  DO 100 i=1,k
2026  m=n+1-i
2027  IF(i.EQ.1) x=a-b/((f+a)*f)
2028  IF(i.EQ.2) x=(z(n)-a)*4d0+z(n)
2029  IF(i.EQ.3) x=(z(n-1)-z(n))*1.6d0+z(n-1)
2030  IF(i.GT.3) x=(z(m+1)-z(m+2))*c+z(m+3)
2031  IF(i.EQ.k.AND.ind.EQ.1) x=0d0
2032  niter=0
2033  check=1d-16
2034  10 pb=1d0
2035  niter=niter+1
2036  IF (niter.LE.100) GO TO 15
2037  check=check*10d0
2038  15 pc=x
2039  dj=a
2040  DO 20 j=2,n
2041  dj=dj+a
2042  pa=pb
2043  pb=pc
2044  20 pc=x*pb+(x*pb-pa)*(dj-a)/dj
2045  pa=a/((pb-x*pc)*f)
2046  pb=pa*pc*(a-x*x)
2047  x=x-pb
2048  IF(dabs(pb).GT.check*dabs(x)) GO TO 10
2049  z(m)=x
2050  w(m)=pa*pa*(a-x*x)
2051  IF(ind1.EQ.0) w(m)=b*w(m)
2052  IF(i.EQ.k.AND.ind.EQ.1) GO TO 100
2053  z(i)=-z(m)
2054  w(i)=w(m)
2055  100 CONTINUE
2056  IF(ind2.NE.1) GO TO 110
2057  print 1100,n
2058  1100 FORMAT(' *** POINTS AND WEIGHTS OF GAUSSIAN QUADRATURE FORMULA',
2059  * ' OF ',i4,'-TH ORDER')
2060  DO 105 i=1,k
2061  zz=-z(i)
2062  105 print 1200,i,zz,i,w(i)
2063  1200 FORMAT(' ',4x,'X(',i4,') = ',f17.14,5x,'W(',i4,') = ',f17.14)
2064  GO TO 115
2065  110 CONTINUE
2066  print 1300,n
2067  1300 FORMAT(' GAUSSIAN QUADRATURE FORMULA OF ',i4,'-TH ORDER IS USED')
2068  115 CONTINUE
2069  IF(ind1.EQ.0) GO TO 140
2070  DO 120 i=1,n
2071  120 z(i)=(a+z(i))/b
2072  140 CONTINUE
2073  RETURN
2074  END
vary
subroutine vary(LAM, MRR, MRI, A, EPS, NP, NGAUSS, X, P, PPI, PIR, PII, R, DR, DDR, DRR, DRI, NMAX)
Definition: ampld.lp.f:978
f1
float f1(float x)
rsp1
subroutine rsp1(X, NG, NGAUSS, REV, EPS, NP, R, DR)
Definition: ampld.lp.f:1030
tmatr0
subroutine tmatr0(NGAUSS, X, W, AN, ANN, S, SS, PPI, PIR, PII, R, DR, DDR, DRR, DRI, NMAX, NCHECK)
Definition: ampld.lp.f:1301
rjb
subroutine rjb(X, Y, U, NMAX, NNMAX)
Definition: ampld.lp.f:1172
vigampl
subroutine vigampl(X, NMAX, M, DV1, DV2)
Definition: ampld.lp.f:860
tmatr
subroutine tmatr(M, NGAUSS, X, W, AN, ANN, S, SS, PPI, PIR, PII, R, DR, DDR, DRR, DRI, NMAX, NCHECK)
Definition: ampld.lp.f:1514
real
#define real
Definition: DbAlgOcean.cpp:26
const
subroutine const(NGAUSS, NMAX, MMAX, P, X, W, AN, ANN, S, SS, NP, EPS)
Definition: ampld.lp.f:924
rsp3
subroutine rsp3(X, NG, NGAUSS, REV, EPS, R, DR)
Definition: ampld.lp.f:1078
rsp4
subroutine rsp4(X, NG, REV, R, DR)
Definition: ampld.lp.f:1112
calcampl
subroutine calcampl(NMAX, LAM, THET0, THET, PHI0, PHI, ALPHA, BETA, S, Z)
Definition: ampld.lp.f:512
vig
subroutine vig(X, NMAX, M, DV1, DV2)
Definition: ampld.lp.f:1793
f2
float f2(float y)
parameter
README for MOD_PR03(V6.1.0) 2. POINTS OF CONTACT it can be either SDP Toolkit or MODIS Packet for Terra input files The orbit validation configuration parameter(LUN 600281) must be either "TRUE" or "FALSE". It needs to be "FALSE" when running in Near Real Time mode
ampl
subroutine ampl(NMAX, DLAM, TL, TL1, PL, PL1, ALPHA, BETA, VV, VH, HV, HH)
Definition: ampld.lp.f:574
sareac
subroutine sareac(EPS, RAT)
Definition: ampld.lp.f:1951
sarea
subroutine sarea(D, RAT)
Definition: ampld.lp.f:1902
surfch
subroutine surfch(N, E, RAT)
Definition: ampld.lp.f:1920
zgetri
subroutine zgetri(N, A, LDA, IPIV, WORK, LWORK, INFO)
Definition: lpd.f:1160
zgetrf
subroutine zgetrf(M, N, A, LDA, IPIV, INFO)
Definition: lpd.f:138
max
#define max(A, B)
Definition: main_biosmap.c:61
tt
subroutine tt(NMAX, NCHECK)
Definition: ampld.lp.f:1852
rsp2
subroutine rsp2(X, NG, REV, EPS, N, R, DR)
Definition: ampld.lp.f:1053
pi
#define pi
Definition: vincenty.c:23
bess
subroutine bess(X, XR, XI, NG, NMAX, NNMAX1, NNMAX2)
Definition: ampld.lp.f:1138
ryb
subroutine ryb(X, Y, V, NMAX)
Definition: ampld.lp.f:1199
f
#define f
Definition: l1_czcs_hdf.c:702
gauss
subroutine gauss(x1, x2, x, w, n)
Definition: 6sm1.f:4055
drop
subroutine drop(RAT)
Definition: ampld.lp.f:1960
calctmat
subroutine calctmat(AXI, RAT, LAM, MRR, MRI, EPS, NP, DDELT, NDGS, NMAX)
Definition: ampld.lp.f:263
cjb
subroutine cjb(XR, XI, YR, YI, UR, UI, NMAX, NNMAX)
Definition: ampld.lp.f:1228
der
subroutine der(T, X, DX)
Definition: der.f:2


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Mission Quality Monitor

HELP

Level 1 & 2 Browser Help

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility Policy

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster