Line data Source code
1 : ! Include shortname defintions, so that the F77 code does not have to be modified to
2 : ! reference the CARMA structure.
3 : #include "carma_globaer.h"
4 :
5 : !! This subroutine computes mie scattering by a stratified sphere,
6 : !! i.e. a particle consisting of a spherical core surrounded by a
7 : !! Spherical shell. The basic code used was that described in the
8 : !! report: " Subroutines for computing the parameters of the
9 : !! electromagnetic radiation scattered by a sphere " J.V. Dave,
10 : !! IBM Scientific Center, Palo Alto , California.
11 : !! Report No. 320 - 3236 .. May 1968 .
12 : !!
13 : !! The modifications for stratified spheres are described in
14 : !! Toon and Ackerman, Appl. Optics, in press, 1981
15 : !!
16 : !! The definitions for the output parameters can be found in "Light
17 : !! scattering by small particles, H.C.Van de Hulst, John Wiley '
18 : !! Sons, Inc., New York, 1957".
19 : !!
20 : !! Also the subroutine computes the capital A function by making use of
21 : !! downward recurrence relationship.
22 : !!
23 : !! @author Brian Toon
24 : !! @version 1981?
25 0 : SUBROUTINE miess(carma,RO,RFR,RFI,THETD,JX,QEXT,QSCAT,QBS,CTBRQS,R,RE2,TMAG2,WVNO,rc)
26 :
27 : ! types
28 : use carma_precision_mod
29 : use carma_constants_mod, only : IT, DEG2RAD
30 : use carma_enums_mod, only : RC_ERROR
31 : use carma_types_mod, only : carma_type
32 : use carma_mod
33 :
34 : implicit none
35 :
36 : type(carma_type), intent(in) :: carma !! the carma object
37 : real(kind=f), intent(in) :: RO !! OUTER (SHELL) RADIUS
38 : real(kind=f), intent(in) :: RFR !! REAL PART OF THE SHELL INDEX OF REFRACTION
39 : real(kind=f), intent(in) :: RFI !! IMAGINARY PART OF THE SHELL INDEX OF REFRACTION
40 : real(kind=f), intent(in) :: R !! CORE RADIUS
41 : real(kind=f), intent(in) :: RE2 !! REAL PART OF THE CORE INDEX OF REFRACTION
42 : real(kind=f), intent(in) :: TMAG2 !! IMAGINARY PART OF THE CORE INDEX OF REFRACTION
43 :
44 : !! ANGLE IN DEGREES BETWEEN THE DIRECTIONS OF THE INCIDENT
45 : !! AND THE SCATTERED RADIATION. THETD(J) IS< OR= 90.0
46 : !! IF THETD(J) SHOULD HAPPEN TO BE GREATER THAN 90.0, ENTER WITH
47 : !! SUPPLEMENTARY VALUE, SEE COMMENTS ON ELTRMX.
48 : real(kind=f), intent(inout) :: THETD(IT)
49 :
50 : !! TOTAL NUMBER OF THETD FOR WHICH THE COMPUTATIONS ARE
51 : !! REQUIRED. JX SHOULD NOT EXCEED IT UNLESS THE DIMENSIONS
52 : !! STATEMENTS ARE APPROPRIATEDLY MODIFIED.
53 : integer, intent(in) :: JX
54 :
55 : real(kind=f), intent(out) :: QEXT !! EFFICIENCY FACTOR FOR EXTINCTION,VAN DE HULST,P.14 ' 127.
56 : real(kind=f), intent(out) :: QSCAT !! EFFICIENCY FACTOR FOR SCATTERING,V.D. HULST,P.14 ' 127.
57 : real(kind=f), intent(out) :: QBS !! BACK SCATTER CROSS SECTION.
58 : real(kind=f), intent(out) :: CTBRQS !! AVERAGE(COSINE THETA) * QSCAT,VAN DE HULST,P.128.
59 : real(kind=f), intent(in) :: WVNO !! 2*PI / WAVELENGTH.
60 : integer, intent(inout) :: rc !! return code, negative indicates failure
61 :
62 : ! Local declarations
63 : real(kind=f), parameter :: EPSILON_MIE = 1.e-14_f
64 :
65 : integer :: I
66 : integer :: J
67 : integer :: K
68 : integer :: M
69 : integer :: N
70 : integer :: NN
71 : integer :: NMX1
72 : integer :: NMX2
73 : integer :: IFLAG
74 : integer :: IACAP
75 :
76 : ! FNAP, FNBP ARE THE PRECEDING VALUES OF FNA, FNB RESPECTIVELY.
77 : complex(kind=f) :: FNAP, FNBP, &
78 : FNA, FNB, RF, RRF, &
79 : RRFX, WM1, FN1, FN2, &
80 : TC1, TC2, WFN(2), Z(4), &
81 : K1, K2, K3, &
82 : RCR, U(8), DH1, &
83 : DH2, DH4, P24H24, P24H21, &
84 : PSTORE, HSTORE, DUMMY, DUMSQ
85 :
86 0 : complex(kind=f), allocatable :: ACAP(:), W(:,:)
87 :
88 :
89 : ! TA(1): REAL PART OF WFN(1). TA(2): IMAGINARY PART OF WFN(1).
90 : ! TA(3): REAL PART OF WFN(2). TA(4): IMAGINARY PART OF WFN(2).
91 : real(kind=f) :: T(5), TA(4), &
92 : PI(3,IT), TAU(3,IT), CSTHT(IT), SI2THT(IT), &
93 : X, X1, X4, Y1, Y4, RX, SINX1, SINX4, COSX1, COSX4, &
94 : EY1, E2Y1, EY4, EY1MY4, EY1PY4, AA, BB, CC, DD, DENOM, &
95 : REALP, AMAGP, QBSR, QBSI, RMM, PIG, RXP4
96 :
97 : !! ELTRMX(I,J,K): ELEMENTS OF THE TRANSFORMATION MATRIX F,V.D.HULST,P.34,45 ' 125.
98 : !! I=1: ELEMENT M SUB 2..I=2: ELEMENT M SUB 1..
99 : !! I = 3: ELEMENT S SUB 21.. I = 4: ELEMENT D SUB 21..
100 : !! ELTRMX(I,J,1) REPRESENTS THE ITH ELEMENT OF THE MATRIX FOR
101 : !! THE ANGLE THETD(J).. ELTRMX(I,J,2) REPRESENTS THE ITH ELEMENT
102 : !! OF THE MATRIX FOR THE ANGLE 180.0 - THETD(J) ..
103 : real(kind=f) :: ELTRMX(4,IT,2)
104 :
105 :
106 : ! IF THE CORE IS SMALL SCATTERING IS COMPUTED FOR THE SHELL ONLY
107 0 : IFLAG = 1
108 0 : if ( R/RO .LT. 1.e-6_f ) IFLAG = 2
109 :
110 0 : if ( JX .gt. IT ) then
111 0 : if (do_print) then
112 0 : write(LUNOPRT, '(a,i3,a)') "miess:: The value of the argument JX=", &
113 0 : JX, " is greater than IT."
114 : end if
115 0 : rc = RC_ERROR
116 0 : return
117 : endif
118 :
119 0 : RF = CMPLX( RFR, -RFI, kind=f )
120 0 : RCR = CMPLX( RE2, -TMAG2, kind=f )
121 0 : X = RO * WVNO
122 0 : K1 = RCR * WVNO
123 0 : K2 = RF * WVNO
124 0 : K3 = CMPLX( WVNO, 0.0_f, kind=f )
125 0 : Z(1) = K2 * RO
126 0 : Z(2) = K3 * RO
127 0 : Z(3) = K1 * R
128 0 : Z(4) = K2 * R
129 0 : X1 = REAL( Z(1) ,f)
130 0 : X4 = REAL( Z(4) ,f)
131 0 : Y1 = aimag( Z(1) )
132 0 : Y4 = aimag( Z(4) )
133 0 : RRF = 1.0_f / RF
134 0 : RX = 1.0_f / X
135 0 : RRFX = RRF * RX
136 0 : T(1) = ( X**2 ) * ( RFR**2 + RFI**2 )
137 0 : T(1) = SQRT( T(1) )
138 0 : NMX1 = 1.10_f * T(1)
139 :
140 : ! The dimension of ACAP.
141 : !
142 : ! In the original program the dimension of ACAP was 7000.
143 : ! For conserving space this should be not much higher than
144 : ! The value, NMX1=1.1*(NREAL**2 + NIMAG**2)**.5 * X + 1
145 0 : IACAP = max(7000, int(1.5_f * NMX1))
146 0 : allocate(ACAP(IACAP))
147 0 : allocate(W(3,IACAP))
148 :
149 0 : NMX2 = T(1)
150 :
151 0 : if ( NMX1 .le. 150 ) then
152 0 : NMX1 = 150
153 0 : NMX2 = 135
154 : endif
155 :
156 0 : ACAP( NMX1+1 ) = ( 0.0_f, 0.0_f )
157 :
158 0 : if ( IFLAG .ne. 2 ) then
159 0 : do N = 1,3
160 0 : W( N,NMX1+1 ) = ( 0.0_f, 0.0_f )
161 : enddo
162 : endif
163 :
164 0 : do N = 1,NMX1
165 0 : NN = NMX1 - N + 1
166 0 : ACAP(NN) = (NN+1) * RRFX - 1.0_f / ( (NN+1) * RRFX + ACAP(NN+1) )
167 0 : if ( IFLAG .ne. 2 ) then
168 0 : do M = 1,3
169 0 : W( M,NN ) = (NN+1) / Z(M+1) - &
170 0 : 1.0_f / ( (NN+1) / Z(M+1) + W( M,NN+1 ) )
171 : enddo
172 : endif
173 : enddo
174 :
175 0 : do J = 1,JX
176 0 : if ( THETD(J) .lt. 0.0_f ) THETD(J) = ABS( THETD(J) )
177 :
178 0 : if ( THETD(J) .le. 0.0_f ) then
179 0 : CSTHT(J) = 1.0_f
180 0 : SI2THT(J) = 0.0_f
181 0 : else if ( THETD(J) .lt. 90.0_f ) then
182 0 : T(1) = THETD(J) * DEG2RAD
183 0 : CSTHT(J) = COS( T(1) )
184 0 : SI2THT(J) = 1.0_f - CSTHT(J)**2
185 0 : else if ( THETD(J) .le. 90.0_f ) then
186 0 : CSTHT(J) = 0.0_f
187 0 : SI2THT(J) = 1.0_f
188 : else
189 0 : if (do_print) then
190 : write(LUNOPRT, '(a,i3)') "miess:: The value of the scattering angle &
191 0 : &is greater than 90.0 Degrees. It is .", THETD(J)
192 : end if
193 0 : rc = RC_ERROR
194 0 : return
195 : end if
196 : enddo
197 :
198 0 : do J = 1,JX
199 0 : PI(1,J) = 0.0_f
200 0 : PI(2,J) = 1.0_f
201 0 : TAU(1,J) = 0.0_f
202 0 : TAU(2,J) = CSTHT(J)
203 : enddo
204 :
205 : ! INITIALIZATION OF HOMOGENEOUS SPHERE
206 0 : T(1) = COS(X)
207 0 : T(2) = SIN(X)
208 0 : WM1 = CMPLX( T(1),-T(2), kind=f )
209 0 : WFN(1) = CMPLX( T(2), T(1), kind=f )
210 0 : TA(1) = T(2)
211 0 : TA(2) = T(1)
212 0 : WFN(2) = RX * WFN(1) - WM1
213 0 : TA(3) = REAL(WFN(2),f)
214 0 : TA(4) = aimag(WFN(2))
215 :
216 0 : if ( IFLAG .ne. 2 ) then
217 0 : N = 1
218 :
219 : ! INITIALIZATION PROCEDURE FOR STRATIFIED SPHERE BEGINS HERE
220 0 : SINX1 = SIN( X1 )
221 0 : SINX4 = SIN( X4 )
222 0 : COSX1 = COS( X1 )
223 0 : COSX4 = COS( X4 )
224 0 : EY1 = EXP( Y1 )
225 0 : E2Y1 = EY1 * EY1
226 0 : EY4 = EXP( Y4 )
227 0 : EY1MY4 = EXP( Y1 - Y4 )
228 0 : EY1PY4 = EY1 * EY4
229 0 : EY1MY4 = EXP( Y1 - Y4 )
230 0 : AA = SINX4 * ( EY1PY4 + EY1MY4 )
231 0 : BB = COSX4 * ( EY1PY4 - EY1MY4 )
232 0 : CC = SINX1 * ( E2Y1 + 1.0_f )
233 0 : DD = COSX1 * ( E2Y1 - 1.0_f )
234 0 : DENOM = 1.0_f + E2Y1 * ( 4.0_f * SINX1 * SINX1 - 2.0_f + E2Y1 )
235 0 : REALP = ( AA * CC + BB * DD ) / DENOM
236 0 : AMAGP = ( BB * CC - AA * DD ) / DENOM
237 0 : DUMMY = CMPLX( REALP, AMAGP, kind=f )
238 0 : AA = SINX4 * SINX4 - 0.5_f
239 0 : BB = COSX4 * SINX4
240 0 : P24H24 = 0.5_f + CMPLX( AA,BB, kind=f ) * EY4 * EY4
241 0 : AA = SINX1 * SINX4 - COSX1 * COSX4
242 0 : BB = SINX1 * COSX4 + COSX1 * SINX4
243 0 : CC = SINX1 * SINX4 + COSX1 * COSX4
244 0 : DD = -SINX1 * COSX4 + COSX1 * SINX4
245 0 : P24H21 = 0.5_f * CMPLX( AA,BB, kind=f ) * EY1 * EY4 + 0.5_f * CMPLX( CC,DD, kind=f ) * EY1MY4
246 0 : DH4 = Z(4) / ( 1.0_f + ( 0.0_f, 1.0_f ) * Z(4) ) - 1.0_f / Z(4)
247 0 : DH1 = Z(1) / ( 1.0_f + ( 0.0_f, 1.0_f ) * Z(1) ) - 1.0_f / Z(1)
248 0 : DH2 = Z(2) / ( 1.0_f + ( 0.0_f, 1.0_f ) * Z(2) ) - 1.0_f / Z(2)
249 0 : PSTORE = ( DH4 + N / Z(4) ) * ( W(3,N) + N / Z(4) )
250 0 : P24H24 = P24H24 / PSTORE
251 0 : HSTORE = ( DH1 + N / Z(1) ) * ( W(3,N) + N / Z(4) )
252 0 : P24H21 = P24H21 / HSTORE
253 0 : PSTORE = ( ACAP(N) + N / Z(1) ) / ( W(3,N) + N / Z(4) )
254 0 : DUMMY = DUMMY * PSTORE
255 0 : DUMSQ = DUMMY * DUMMY
256 :
257 : ! NOTE: THE DEFINITIONS OF U(I) IN THIS PROGRAM ARE NOT THE SAME AS
258 : ! THE USUBI DEFINED IN THE ARTICLE BY TOON AND ACKERMAN. THE
259 : ! CORRESPONDING TERMS ARE:
260 : ! USUB1 = U(1) USUB2 = U(5)
261 : ! USUB3 = U(7) USUB4 = DUMSQ
262 : ! USUB5 = U(2) USUB6 = U(3)
263 : ! USUB7 = U(6) USUB8 = U(4)
264 : ! RATIO OF SPHERICAL BESSEL FTN TO SPHERICAL HENKAL FTN = U(8)
265 :
266 0 : U(1) = K3 * ACAP(N) - K2 * W(1,N)
267 0 : U(2) = K3 * ACAP(N) - K2 * DH2
268 0 : U(3) = K2 * ACAP(N) - K3 * W(1,N)
269 0 : U(4) = K2 * ACAP(N) - K3 * DH2
270 0 : U(5) = K1 * W(3,N) - K2 * W(2,N)
271 0 : U(6) = K2 * W(3,N) - K1 * W(2,N)
272 0 : U(7) = ( 0.0_f, -1.0_f ) * ( DUMMY * P24H21 - P24H24 )
273 0 : U(8) = TA(3) / WFN(2)
274 :
275 : FNA = U(8) * ( U(1)*U(5)*U(7) + K1*U(1) - DUMSQ*K3*U(5) ) / &
276 0 : ( U(2)*U(5)*U(7) + K1*U(2) - DUMSQ*K3*U(5) )
277 : FNB = U(8) * ( U(3)*U(6)*U(7) + K2*U(3) - DUMSQ*K2*U(6) ) / &
278 0 : ( U(4)*U(6)*U(7) + K2*U(4) - DUMSQ*K2*U(6) )
279 : else
280 0 : TC1 = ACAP(1) * RRF + RX
281 0 : TC2 = ACAP(1) * RF + RX
282 0 : FNA = ( TC1 * TA(3) - TA(1) ) / ( TC1 * WFN(2) - WFN(1) )
283 0 : FNB = ( TC2 * TA(3) - TA(1) ) / ( TC2 * WFN(2) - WFN(1) )
284 : endif
285 :
286 0 : FNAP = FNA
287 0 : FNBP = FNB
288 0 : T(1) = 1.50_f
289 :
290 : ! FROM HERE TO THE STATMENT NUMBER 90, ELTRMX(I,J,K) HAS
291 : ! FOLLOWING MEANING:
292 : ! ELTRMX(1,J,K): REAL PART OF THE FIRST COMPLEX AMPLITUDE.
293 : ! ELTRMX(2,J,K): IMAGINARY PART OF THE FIRST COMPLEX AMPLITUDE.
294 : ! ELTRMX(3,J,K): REAL PART OF THE SECOND COMPLEX AMPLITUDE.
295 : ! ELTRMX(4,J,K): IMAGINARY PART OF THE SECOND COMPLEX AMPLITUDE.
296 : ! K = 1 : FOR THETD(J) AND K = 2 : FOR 180.0 - THETD(J)
297 : ! DEFINITION OF THE COMPLEX AMPLITUDE: VAN DE HULST,P.125.
298 0 : FNA = T(1) * FNA
299 0 : FNB = T(1) * FNB
300 :
301 0 : do J = 1,JX
302 0 : TC1 = FNA * PI(2,J) + FNB * TAU(2,J)
303 0 : TC2 = FNB * PI(2,J) + FNA * TAU(2,J)
304 0 : ELTRMX(1,J,1) = real(TC1)
305 0 : ELTRMX(2,J,1) = aimag(TC1)
306 0 : ELTRMX(3,J,1) = real(TC2)
307 0 : ELTRMX(4,J,1) = aimag(TC2)
308 0 : TC1 = FNA * PI(2,J) - FNB * TAU(2,J)
309 0 : TC2 = FNB * PI(2,J) - FNA * TAU(2,J)
310 0 : ELTRMX(1,J,2) = real(TC1)
311 0 : ELTRMX(2,J,2) = aimag(TC1)
312 0 : ELTRMX(3,J,2) = real(TC2)
313 0 : ELTRMX(4,J,2) = aimag(TC2)
314 : enddo
315 :
316 0 : QEXT = 2.0_f * ( real(FNA) + real(FNB) )
317 : QSCAT = ( real(FNA)**2 + aimag(FNA)**2 + &
318 0 : real(FNB)**2 + aimag(FNB)**2 ) / 0.75_f
319 0 : CTBRQS = 0.0_f
320 0 : TC1 = -2.0_f * (FNB - FNA)
321 0 : QBSR = real(TC1)
322 0 : QBSI = aimag(TC1)
323 0 : RMM = -1.0_f
324 0 : N = 2
325 :
326 : ! Iterate until the answer converges.
327 0 : T(4) = EPSILON_MIE
328 :
329 0 : do while ( T(4) .ge. EPSILON_MIE )
330 :
331 0 : T(1) = 2*N - 1
332 0 : T(2) = N - 1
333 0 : T(3) = 2*N + 1
334 :
335 0 : do J = 1,JX
336 0 : PI(3,J) = ( T(1) * PI(2,J) * CSTHT(J) - N * PI(1,J) ) / T(2)
337 : TAU(3,J) = CSTHT(J) * ( PI(3,J) - PI(1,J) ) - &
338 0 : T(1) * SI2THT(J) * PI(2,J) + TAU(1,J)
339 : end do
340 :
341 : ! HERE SET UP HOMOGENEOUS SPHERE
342 0 : WM1 = WFN(1)
343 0 : WFN(1) = WFN(2)
344 0 : TA(1) = REAL(WFN(1),f)
345 0 : TA(2) = aimag(WFN(1))
346 0 : TA(4) = aimag(WFN(2))
347 0 : WFN(2) = T(1) * RX * WFN(1) - WM1
348 0 : TA(3) = REAL(WFN(2),f)
349 :
350 0 : if ( IFLAG .ne. 2 ) then
351 :
352 : ! HERE SET UP STRATIFIED SPHERE
353 0 : DH2 = - N / Z(2) + 1.0_f / ( N / Z(2) - DH2 )
354 0 : DH4 = - N / Z(4) + 1.0_f / ( N / Z(4) - DH4 )
355 0 : DH1 = - N / Z(1) + 1.0_f / ( N / Z(1) - DH1 )
356 0 : PSTORE = ( DH4 + N / Z(4) ) * ( W(3,N) + N / Z(4) )
357 0 : P24H24 = P24H24 / PSTORE
358 0 : HSTORE = ( DH1 + N / Z(1) ) * ( W(3,N) + N / Z(4) )
359 0 : P24H21 = P24H21 / HSTORE
360 0 : PSTORE = ( ACAP(N) + N / Z(1) ) / ( W(3,N) + N / Z(4) )
361 0 : DUMMY = DUMMY * PSTORE
362 0 : DUMSQ = DUMMY * DUMMY
363 :
364 0 : U(1) = K3 * ACAP(N) - K2 * W(1,N)
365 0 : U(2) = K3 * ACAP(N) - K2 * DH2
366 0 : U(3) = K2 * ACAP(N) - K3 * W(1,N)
367 0 : U(4) = K2 * ACAP(N) - K3 * DH2
368 0 : U(5) = K1 * W(3,N) - K2 * W(2,N)
369 0 : U(6) = K2 * W(3,N) - K1 * W(2,N)
370 0 : U(7) = ( 0.0_f, -1.0_f ) * ( DUMMY * P24H21 - P24H24 )
371 0 : U(8) = TA(3) / WFN(2)
372 :
373 : FNA = U(8) * ( U(1)*U(5)*U(7) + K1*U(1) - DUMSQ*K3*U(5) ) / &
374 0 : ( U(2)*U(5)*U(7) + K1*U(2) - DUMSQ*K3*U(5) )
375 : FNB = U(8) * ( U(3)*U(6)*U(7) + K2*U(3) - DUMSQ*K2*U(6) ) / &
376 0 : ( U(4)*U(6)*U(7) + K2*U(4) - DUMSQ*K2*U(6) )
377 : endif
378 :
379 0 : TC1 = ACAP(N) * RRF + N * RX
380 0 : TC2 = ACAP(N) * RF + N * RX
381 0 : FN1 = ( TC1 * TA(3) - TA(1) ) / ( TC1 * WFN(2) - WFN(1) )
382 0 : FN2 = ( TC2 * TA(3) - TA(1) ) / ( TC2 * WFN(2) - WFN(1) )
383 0 : M = WVNO * R
384 :
385 0 : if ( N .ge. M ) then
386 0 : if ( IFLAG .ne. 2 ) then
387 0 : if ( abs( ( FN1-FNA ) / FN1 ) .LT. EPSILON_MIE .AND. &
388 : abs( ( FN2-FNB ) / FN2 ) .LT. EPSILON_MIE ) IFLAG = 2
389 :
390 0 : if ( IFLAG .ne. 1 ) then
391 : FNA = FN1
392 : FNB = FN2
393 : endif
394 : else
395 : FNA = FN1
396 : FNB = FN2
397 : endif
398 : endif
399 :
400 0 : T(5) = N
401 0 : T(4) = T(1) / ( T(5) * T(2) )
402 0 : T(2) = ( T(2) * ( T(5) + 1.0_f ) ) / T(5)
403 :
404 : CTBRQS = CTBRQS + T(2) * ( real(FNAP) * real(FNA) + &
405 : aimag(FNAP) *aimag(FNA) &
406 : + real(FNBP) * real(FNB) + &
407 : aimag(FNBP) *aimag(FNB) ) &
408 : + T(4) * ( real(FNAP) * real(FNBP) + &
409 0 : aimag(FNAP) *aimag(FNBP) )
410 0 : QEXT = QEXT + T(3) * ( real(FNA) + real(FNB) )
411 :
412 : ! $ T(3), real(FNA), real(FNB), QEXT
413 : T(4) = real(FNA)**2 + aimag(FNA)**2 + &
414 0 : real(FNB)**2 + aimag(FNB)**2
415 0 : QSCAT = QSCAT + T(3) * T(4)
416 0 : RMM = -RMM
417 0 : TC1 = T(3)*RMM*(FNB - FNA)
418 0 : QBSR = QBSR + real(TC1)
419 0 : QBSI = QBSI + aimag(TC1)
420 :
421 0 : T(2) = N * (N+1)
422 0 : T(1) = T(3) / T(2)
423 0 : K = (N/2)*2
424 :
425 0 : do J = 1,JX
426 0 : TC1 = FNA * PI(3,J) + FNB * TAU(3,J)
427 0 : TC2 = FNB * PI(3,J) + FNA * TAU(3,J)
428 0 : ELTRMX(1,J,1) = ELTRMX(1,J,1)+T(1)* real(TC1)
429 0 : ELTRMX(2,J,1) = ELTRMX(2,J,1)+T(1)*aimag(TC1)
430 0 : ELTRMX(3,J,1) = ELTRMX(3,J,1)+T(1)* real(TC2)
431 0 : ELTRMX(4,J,1) = ELTRMX(4,J,1)+T(1)*aimag(TC2)
432 :
433 0 : IF ( K .EQ. N ) THEN
434 0 : TC1 = -FNA * PI(3,J) + FNB * TAU(3,J)
435 0 : TC2 = -FNB * PI(3,J) + FNA * TAU(3,J)
436 : ELSE
437 0 : TC1 = FNA * PI(3,J) - FNB * TAU(3,J)
438 0 : TC2 = FNB * PI(3,J) - FNA * TAU(3,J)
439 : END IF
440 0 : ELTRMX(1,J,2) = ELTRMX(1,J,2)+T(1)* real(TC1)
441 0 : ELTRMX(2,J,2) = ELTRMX(2,J,2)+T(1)*aimag(TC1)
442 0 : ELTRMX(3,J,2) = ELTRMX(3,J,2)+T(1)* real(TC2)
443 0 : ELTRMX(4,J,2) = ELTRMX(4,J,2)+T(1)*aimag(TC2)
444 :
445 : enddo
446 :
447 0 : if ( T(4) .ge. EPSILON_MIE ) then
448 0 : N = N + 1
449 :
450 0 : do J = 1,JX
451 0 : PI(1,J) = PI(2,J)
452 0 : PI(2,J) = PI(3,J)
453 0 : TAU(1,J) = TAU(2,J)
454 0 : TAU(2,J) = TAU(3,J)
455 : enddo
456 :
457 0 : FNAP = FNA
458 0 : FNBP = FNB
459 :
460 0 : if ( N .gt. NMX2 ) then
461 0 : if (do_print) write(LUNOPRT, '(a)') "miess:: The upper limit for acap is not enough."
462 0 : rc = RC_ERROR
463 0 : return
464 : endif
465 : endif
466 : enddo
467 :
468 : ! Calculate the results.
469 0 : do J = 1,JX
470 0 : do K = 1,2
471 0 : do I= 1,4
472 0 : T(I) = ELTRMX(I,J,K)
473 : enddo
474 :
475 0 : ELTRMX(2,J,K) = T(1)**2 + T(2)**2
476 0 : ELTRMX(1,J,K) = T(3)**2 + T(4)**2
477 0 : ELTRMX(3,J,K) = T(1) * T(3) + T(2) * T(4)
478 0 : ELTRMX(4,J,K) = T(2) * T(3) - T(4) * T(1)
479 : enddo
480 : enddo
481 :
482 0 : T(1) = 2.0_f * RX**2
483 0 : QEXT = QEXT * T(1)
484 0 : QSCAT = QSCAT * T(1)
485 0 : CTBRQS = 2.0_f * CTBRQS * T(1)
486 :
487 : ! QBS IS THE BACK SCATTER CROSS SECTION
488 0 : PIG = ACOS(-1.0_f)
489 0 : RXP4 = RX*RX/(4.0_f*PIG)
490 0 : QBS = RXP4*(QBSR**2 + QBSI**2)
491 :
492 0 : deallocate(ACAP)
493 0 : deallocate(W)
494 :
495 0 : return
496 0 : end
|
