Line data Source code
1 : !! This module (fractal_meanfield_mod.F90) contains the main routines
2 : !! necessary to calculate the solution of the mean field approximation
3 : !! for a dry fractal particle composed of identical spherical monomers.
4 : !! This is used to generate optical properties for these paticles in CARMA.
5 : !!
6 : !! See Botet et al. 1997 "Mean-field approximation of Mie
7 : !! scattering by fractal aggregates of identical spheres."
8 : !! Applied Optics 36(33) 8791-8797
9 : !!
10 : !! Original code from P. Rannou and R. Botet.
11 : !! Translated to F90 and ported into CARMA by E. Wolf
12 : !!
13 : !! master: fractal_meanfield calling: cmie,ludcmpc,lubksbc,dqagi
14 : !!
15 : !! calculating the monomer Mie scattering
16 : !! - SUBROUTINE cmie() calling: intmie()
17 : !! - SUBROUTINE intmie() calling: intmie()
18 : !!
19 : !! calculating the matrix elements
20 : !! - FUNCTION funa() calling: dqag,fpl
21 : !! - FUNCTION fpl() calling: calcpoly
22 : !! - FUNCTION calcpoly()
23 : !! - FUNCTION funb_n()
24 : !! - FUNCTION funs_n() calling: dq2agi,xfreal_n,xfimag_n
25 : !! - FUNCTION xfreal_n() calling: besseljy,phi
26 : !! - FUNCTION xfimag_n() calling: besseljy,phi
27 : !! - FUNCTION BESSELJY()
28 : !!
29 : !! Routines to calculate the scattered wave
30 : !! of monomer:
31 : !! - FUNCTION fpi() calling: calcpoly()
32 : !! - FUNCTION ftau() calling: calcpoly()
33 : !! of agglomerate/cluster:
34 : !! - FUNCTION fp1()
35 : !!
36 : !! Routines related to the probability distribution:
37 : !! - FUNCTION anorm() calling: dqdagi,fdval
38 : !! - FUNCTION fdval()
39 : !! - FUNCTION phi() calling: fdval
40 : !! - FUNCTION fco() calling: fdval
41 : !!
42 : !! @author P. Rannou, R. Botet, Eric Wolf
43 : !! version March 2013
44 : module fractal_meanfield_mod
45 :
46 : use carma_precision_mod
47 : use carma_enums_mod
48 : use carma_constants_mod
49 : use carma_types_mod
50 : use carma_mod
51 :
52 : use adgaquad_types_mod
53 : use adgaquad_mod
54 :
55 : implicit none
56 :
57 : private
58 :
59 : public :: fractal_meanfield
60 :
61 : ! Private module varibles: Moved from COMMON blocks
62 : integer, parameter :: nmi=40
63 : integer, parameter :: n2m = 2*nmi
64 :
65 : contains
66 :
67 : !!
68 : !! Generate optical properties for CARMA fractal particles.
69 : !!
70 : !! See Botet et al. 1997 "Mean-field approximation of Mie
71 : !! scattering by fractal aggregates of identical spheres."
72 : !! Applied Optics 36(33) 8791-8797
73 : !!
74 : !! @author P.Rannou, R.Botet, Eric Wolf
75 : !! @version March 2013
76 0 : subroutine fractal_meanfield(carma, xl_in, xk_in, xn_in, nb_in, alpha_in, &
77 : df_in, rmon,xv, ang, Qext, Qsca, gfac, rc)
78 :
79 : ! some of these may be included in carma object
80 : type(carma_type), intent(in) :: carma !! the carma object
81 : real(kind=f),intent(in) :: xl_in !! Wavelength [microns]
82 : real(kind=f),intent(in) :: xk_in !! imaginary index of refraction
83 : real(kind=f),intent(in) :: xn_in !! real index of refraction
84 : real(kind=f),intent(in) :: nb_in !! number of monomers
85 : real(kind=f),intent(in) :: alpha_in !! Packing coefficient
86 : real(kind=f),intent(in) :: df_in !! Fractal dimension
87 : real(kind=f),intent(in) :: rmon !! monomer size [microns]
88 : real(kind=f),intent(in) :: xv !! set to 1
89 : real(kind=f),intent(in) :: ang !! angle set to zero
90 : real(kind=f),intent(out) :: Qext !! EFFICIENCY FACTOR FOR EXTINCTION
91 : real(kind=f),intent(out) :: Qsca !! EFFICIENCY FACTOR FOR SCATTERING
92 : real(kind=f),intent(out) :: gfac !! asymmetry factor
93 : integer,intent(inout) :: rc !! return code, negative indicates failure
94 :
95 : ! Local declarations
96 : integer, parameter :: nth = 10001
97 : integer, parameter :: maxsub = 50000
98 : integer, parameter :: lenw = 200000
99 : integer :: nstop ! index with the last mie-coefficient
100 : integer :: n1stop
101 : integer :: pp,tt,mm
102 : real(kind=f) :: krg,rg ! Particle structure
103 : real(kind=f) :: sigmas,sigmae,nc1
104 : real(kind=f) :: sigext ! extinction cross section
105 : real(kind=f) :: sigsca ! scattering cross section
106 : real(kind=f) :: sigabs ! absorption cross section
107 : real(kind=f) :: totg ! asymmetry parameter
108 : real(kind=f) :: sigext2,sigext3
109 : real(kind=f) :: rems ! radius of equivalent mass sphere
110 : real(kind=f) :: gems ! geometric cross-section of equivalent mass sphere
111 : real(kind=f) :: dthetar,angler,weight
112 : real(kind=f) :: sumsca
113 : real(kind=f) :: lbd,beta ! optical characteristics
114 : real(kind=f) :: sstest(0:nf)
115 : real(kind=f) :: setest(0:nf)
116 : real(kind=f) :: xl(39) ! place holder for wavelength
117 : real(kind=f) :: xn(39) ! place holder for real index of refraction
118 : real(kind=f) :: xk(39) ! place holder for imaginary index of refraction
119 : real(kind=f) :: val
120 : real(kind=f) :: funca(nmi,nmi,0:n2m) ! for storage of funa(nu,n;p)
121 : complex(kind=f) :: res ! for storage of funs_n
122 : complex(kind=f) :: funcs(0:n2m)
123 : real(kind=f) :: s11(0:nth-1)
124 : real(kind=f) :: s11_n(0:nth-1)
125 : real(kind=f) :: xint(0:nth-1)
126 : real(kind=f) :: wom
127 : real(kind=f) :: pol(0:nth-1)
128 : complex(kind=f) :: s01,s02
129 : complex(kind=f) :: s1(0:nth-1)
130 : complex(kind=f) :: s2(0:nth-1)
131 : complex(kind=f) :: ajt
132 : complex(kind=f) :: an(nf)
133 : complex(kind=f) :: bn(nf)
134 : complex(kind=f) :: ni,i,id,onec,zeroc
135 : complex(kind=f) :: d1(nmi)
136 : complex(kind=f) :: d2(nmi)
137 : complex(kind=f) :: Ap1(nmi,nmi)
138 : complex(kind=f) :: Bp1(nmi,nmi)
139 : complex(kind=f) :: cvec(n2m)
140 : complex(kind=f) :: EpABC(n2m,n2m)
141 : integer :: luindx(n2m) ! For LU decomposition
142 : real(kind=f) :: dlu
143 : integer :: ifail
144 : integer :: iwork(maxsub)
145 : integer :: neval,nsubin
146 : real(kind=f) :: work(lenw)
147 :
148 : ! Previously these were implicitly defined
149 : real(kind=f) :: angle, pi, rn, ri
150 : real(kind=f) :: deltas, deltae, xfact
151 : real(kind=f) :: bound, errrel, p1, dp1
152 : real(kind=f) :: errabs, total
153 : integer :: n2stop, n3stop, ntheta, ii, kk, nn, jj, iy, ir, q, interv
154 : real(kind=f) :: a0, c0, a1, c1, a2, c2
155 : real(kind=f) :: qabs
156 :
157 : ! Previously these were globals, which wouldn't be thread safe.
158 : type(adgaquad_vars_type) :: fx_vars
159 :
160 : integer :: info
161 :
162 : external ZGESV ! lapack routine to solve complex matrix eqn.
163 :
164 0 : errabs=0.0_f
165 :
166 : ! Set the return code to default to okay.
167 0 : rc = RC_OK
168 :
169 : ! *** Set from input arguments
170 0 : fx_vars%nb = nb_in
171 0 : fx_vars%df = df_in
172 0 : fx_vars%alpha = alpha_in
173 0 : xl(1) = xl_in
174 0 : xk(1) = xk_in
175 0 : xn(1) = xn_in
176 :
177 : ! *** Complex constants 1, 1, identity(1,1), zero(0,0) :
178 0 : i = cmplx(0._f,1._f,kind=f)
179 0 : onec = cmplx(1._f,0._f,kind=f)
180 0 : id = cmplx(1._f,1._f,kind=f)
181 0 : zeroc = cmplx(0._f,0._f,kind=f)
182 :
183 : ! Other initializations
184 0 : funca(:,:,:) = 0.0_f
185 0 : fx_vars%a = rmon *1.e-2_f ! a = r_monomer in m
186 0 : beta=ang*(3.1415926_f / 180._f) ! =0 when ang=0
187 0 : Ap1(:,:) = zeroc
188 0 : Bp1(:,:) = zeroc
189 0 : sstest(:) = 0.0_f
190 0 : setest(:) = 0.0_f
191 :
192 : ! ****************************************************************
193 : ! *** Definition and calculation of factorials 0 - nf
194 : ! *** (nf set in adgaqaud_types_mod.F90)
195 : ! *** and storage [ real*8 fact() (double prec.) ]
196 : ! ****************************************************************
197 :
198 0 : fx_vars%fact(0)=1._f ! factorials fact(n)=n!
199 0 : do ii=1,nf
200 0 : fx_vars%fact(ii) = fx_vars%fact(ii-1)*ii*1._f
201 : end do
202 :
203 0 : pi=4._f*atan(1._f) ! 3.1415926535
204 0 : fx_vars%coeff=anorm(carma,fx_vars,rc)
205 0 : if (rc < 0) return
206 :
207 : ! ****************************************************************
208 : ! anorm() integrated INT_0^inf[ x**(df-1.)*exp(-x**df/2._f) dx ]
209 : ! and occupied
210 : ! anorm := 4 pi * INT_0^inf[ x**(df-1.)*exp(-x**df/2._f) dx ]
211 : ! == geometric scalingfactor Eq.(10) in [Botet et al, 1995]
212 : ! c := 0.5
213 : ! ****************************************************************
214 :
215 0 : kk=1
216 0 : ni=xn(kk)*1._f+i*xk(kk)*xv*1._f ! ni := complex index of refraction of monomer
217 : ! (xv := 1 ; input parameter in file "calpha")
218 0 : lbd=xl(kk)*1.e-6_f ! lbd := wavelength in m
219 : ! (in matrix medium / material !)
220 0 : fx_vars%k=2._f*pi/lbd ! k := abs.val. of wavevector in m^-1
221 : ! (in matrix medium / material !)
222 :
223 : ! *** ******************************************************************
224 : ! *** Calculation of Mie coefficients for monomer scattering
225 : ! *** up to a maximum order of nf=50
226 : ! *** ******************************************************************
227 :
228 0 : do ii=1,nf
229 0 : an(ii) = zeroc
230 0 : bn(ii) = zeroc
231 : end do
232 :
233 0 : rn=xn(kk) ! Re(relative_n_complex,monomer)
234 0 : ri=xk(kk)*xv ! Im(relative_n_complex,monomer)
235 : ! xv should be set to 1 (sse above)
236 : ! a = monomer sphere radius
237 : ! lbd = wavelength in matrix medium
238 :
239 : ! Call Mie routine
240 0 : call cmie(lbd,rn,ri,fx_vars%a,an,bn,nstop)
241 :
242 0 : do ii=1,nf
243 0 : if (an(ii).ne.0._f) nstop=ii
244 : end do
245 :
246 : ! nstop is now the index with the last mie-coefficient
247 : ! (highest index i) an(i) not equal zero.
248 : ! since all the an were set to zero before calling
249 : ! cmie(), nstop is the termination index used in cmie()
250 : ! or in intmie(). Usually, a termination index
251 : ! nstop = INT( 2 + x + 4 x^(1/3) ) is used; in intmie(),
252 : ! however, a value of
253 : ! nstop := MAX( INT(...), |m*x| )+15 is used !???
254 :
255 0 : sigmas=0._f
256 0 : sigmae=0._f
257 :
258 0 : do nn=1,nstop
259 0 : nc1=abs(an(nn))**2._f+abs(bn(nn))**2._f
260 0 : nc1=nc1*(2._f*nn+1._f)
261 0 : sigmas=sigmas+nc1*(2._f*3.14159265_f)/(fx_vars%k**2._f)
262 0 : nc1=real(an(nn)+bn(nn))
263 0 : nc1=nc1*(2._f*nn+1._f)
264 0 : sigmae=sigmae+nc1*(2._f*3.14159265_f)/(fx_vars%k**2._f)
265 0 : sstest(nn)=sigmas
266 0 : setest(nn)=sigmae
267 0 : deltas=abs(sstest(nn-1)-sstest(nn))/sstest(nn)
268 0 : deltae=abs(setest(nn-1)-setest(nn))/setest(nn)
269 0 : if(deltas.gt.1.e-6_f) n2stop=nn
270 0 : if(deltae.gt.1.e-6_f) n3stop=nn
271 : end do
272 :
273 0 : n1stop=n2stop
274 0 : if (n3stop.gt.n2stop) n1stop=n3stop
275 : ! The order of the set of linear equations is chosen
276 : ! as the number of mie coefficients where the sum yielding
277 : ! the monomer ext./scatt. cross sections do not change more
278 : ! than 1.D-3 compared to the values with one summand less.
279 :
280 0 : rg=fx_vars%alpha*fx_vars%nb**(1._f/fx_vars%df)*fx_vars%a ! rg := radius of gyration
281 0 : krg=fx_vars%k*rg
282 : ntheta=7800 !180-int(krg**.5*28*log10(dxk(kk)))
283 : if (ntheta*0.5_f .eq. (ntheta/2)*1._f) ntheta=ntheta+1
284 :
285 : ! *** ******************************************************************
286 : ! *** FIRST PART: Solution of self consistent mean field equation,
287 : ! *** i.e. the set of linear equations (SLE) defining the
288 : ! *** mean field coefficients d^1_1,n and d^2_1,n
289 : ! *** according to Eq.(12) of (Botet 1997)
290 : ! *** To do so,
291 : ! *** - matrix elements A^1,nu_1,n and B^1,nu_1,n
292 : ! *** are calculated with Eq.(13), using Eqns.(14)-(16)
293 : ! *** - the set of lin. Eqns. is solved yielding the d's
294 : ! *** ******************************************************************
295 : ! *** Eq.(12) of Botet et al 1997 defines a matrix eqn. of order 2N :
296 : ! *** (since N=n1stop, 2N = 2 * n1stop = order of SLE)
297 : ! ***
298 : ! *** EpABC * dvec = cvec with
299 : ! ***
300 : ! *** dvec and cvec being the 2N-vectors
301 : ! ***
302 : ! *** ( d^(1)_1,1 ) ( a_1 )
303 : ! *** ( d^(1)_1,2 ) ( a_2 )
304 : ! *** ( ... ) ( ... )
305 : ! *** ( d^(1)_1,n ) ( a_n )
306 : ! *** dvec := ( d^(2)_1,1 ) and cvec := ( b_1 ) and further
307 : ! *** ( d^(2)_1,2 ) ( b_2 )
308 : ! *** ( ... ) ( ... )
309 : ! *** ( d^(2)_1,n ) ( b_n )
310 : ! ***
311 : ! ***
312 : ! *** EpABC := 1 + AB * C where AB, 1, C are the 2N*2N - matrices
313 : ! ***
314 : ! *** ( a_1 0 0 ... ... 0 )
315 : ! *** ( 0 a_2 0 ... ... 0 )
316 : ! *** AB := ( 0 0 ... 0 0 0 ... 0 ) 1 := 2N*2N unity matrix
317 : ! *** ( 0 0 ... a_n 0 0 ... 0 )
318 : ! *** ( 0 0 ... 0 b_1 0 ... 0 )
319 : ! *** ( ... ... ... 0 )
320 : ! *** ( ... ... 0 b_n-1 0 )
321 : ! *** ( 0 0 ... ... 0 0 b_n)
322 : ! ***
323 : ! *** and ( A B )
324 : ! *** C := ( ) where A and B are the two N*N matrices
325 : ! *** ( B A ) given by Eq.(13),
326 : ! *** including the factor (N_monomers - 1):
327 : ! ***
328 : ! *** A_n,nu := (N_m-1) * A_(1,n)^(1,nu) and
329 : ! *** B_n,nu := (N_m-1) * B_(1,n)^(1,nu)
330 : ! ***
331 : ! *** (A_(1,n... and B_(1,n... according to Eq.(13) of Botet 1997)
332 : ! *** ******************************************************************
333 :
334 0 : n2stop = 2 * n1stop ! n2stop = order of SLE
335 :
336 : ! *** ******************************************************************
337 : ! *** Error handling moved from xfreal_n, xfimag_n. Calculations fail
338 : ! *** in integration package when n2stop > 48. n2stop is related to the
339 : ! *** number of complex mie scattering coefficients used in teh calculation
340 : ! *** which is in turn related to the size parameter of monomers.
341 : ! *** If nstop>48 end calculation here instead of continuing.
342 :
343 0 : if (n2stop.gt.48) then
344 0 : if (carma%f_do_print) then
345 : write(carma%f_LUNOPRT, *) "fractal_meanfield_mod::n2stop greater &
346 0 : &than 48. Size parameter (2*pi*rmon/lambda): ", &
347 0 : 2._f*3.14159265_f*fx_vars%a/lbd, "Monomer Size parameter &
348 0 : &must be less than ~17."
349 : end if
350 0 : rc = RC_ERROR
351 0 : return
352 : endif
353 :
354 :
355 : ! *** ******************************************************************
356 0 : do ii=1,n1stop
357 0 : cvec(ii) = an(ii) ! right hand side vector
358 0 : cvec(n1stop+ii) = bn(ii)
359 : end do
360 :
361 0 : do pp=0,n2stop !variable p
362 0 : res = funs_n(carma,fx_vars,pp,rc) ! Eq.(16) S_p(k R_g)
363 0 : if (rc < 0) return
364 0 : funcs(pp) = res
365 : end do
366 :
367 : ! Calculate terms A and B
368 :
369 : ! *** loops over indices nu,n,p :
370 0 : do ii=1,n1stop !variable n
371 0 : do jj=1,n1stop !variable nu
372 0 : mm = IABS(ii-jj)
373 0 : tt = ii+jj
374 :
375 : ! *** ******************************************************************
376 : ! calculation of A_(1,n)^(1,nu) according to Eq.(13)
377 : ! *** ******************************************************************
378 0 : do pp=mm,tt
379 :
380 0 : funca(jj,ii,pp) = funa(carma,fx_vars,jj,ii,pp,rc) ! Eq.(14) a(nu,n;p)a
381 0 : if (rc < 0) return
382 :
383 0 : Ap1(ii,jj) = Ap1(ii,jj) + &
384 : ( onec * (ii*(ii+1)+jj*(jj+1)-pp*(pp+1)) ) &
385 0 : * funca(jj,ii,pp) * funcs(pp)
386 : end do ! loop over pp (variable p)
387 :
388 : ! scaling factors of eq.(13), factor (N_mon-1) from eq.(12)
389 0 : Ap1(ii,jj) = Ap1(ii,jj) * (2._f*jj+1._f)/(jj*(jj*1._f+1._f))
390 0 : Ap1(ii,jj) = Ap1(ii,jj) * (fx_vars%nb-1._f) / (ii*(ii*1._f+1._f))
391 :
392 : ! *** ******************************************************************
393 : ! calculation of B_(1,n)^(1,nu) according to Eq.(13)
394 : ! *** ******************************************************************
395 0 : do pp=mm,tt
396 0 : Bp1(ii,jj) = Bp1(ii,jj) + funb_n(jj,ii,pp,funca) * funcs(pp)
397 : end do ! loop over pp (variable p)
398 :
399 : ! scaling factors of eq.(13), factor (N_mon-1) from eq.(12)
400 0 : Bp1(ii,jj) = Bp1(ii,jj) * (2._f*jj+1._f)/(jj*(jj*1._f+1._f))
401 0 : Bp1(ii,jj) = Bp1(ii,jj) * (fx_vars%nb-1._f) * 2._f/(ii*(ii*1._f+1._f))
402 : end do ! loop over jj=1,n1stop (variable nu)
403 : end do ! loop over ii=1,n1stop (variable n)
404 :
405 : ! *** ******************************************************************
406 : ! End of Calculation of terms A and B
407 :
408 : ! *** ******************************************************************
409 : ! *** Setup and solution of matrix equation of order 2N ( = n2stop )
410 : ! *** constituted by eq.(12)
411 : ! *** ******************************************************************
412 : ! *** matrix product (AB * C) (definitions see above)
413 0 : do ii=1,n1stop
414 0 : do jj=1,n1stop
415 0 : EpABC(ii,jj) = an(ii) * Ap1(ii,jj)
416 0 : EpABC(ii,jj+n1stop) = an(ii) * Bp1(ii,jj)
417 0 : EpABC(ii+n1stop,jj) = bn(ii) * Bp1(ii,jj)
418 0 : EpABC(ii+n1stop,jj+n1stop) = bn(ii) * Ap1(ii,jj)
419 : end do
420 : end do
421 :
422 : ! *** ******************************************************************
423 : ! *** add 2N*2N unity matrix
424 0 : do ii=1,n1stop
425 0 : EpABC(ii,ii) = EpABC(ii,ii) + onec
426 0 : EpABC(ii+n1stop,ii+n1stop) = EpABC(ii+n1stop,ii+n1stop) + onec
427 : end do
428 :
429 : ! ======================================================================
430 : ! *** solve matrix equation using lapack routine
431 0 : call ZGESV(n2stop, 1, EpABC, n2m, luindx, cvec, n2m, info)
432 0 : if (info/=0) then
433 0 : rc = info
434 0 : return
435 : end if
436 :
437 0 : do ii=1,n1stop
438 0 : d1(ii) = cvec(ii)
439 0 : d2(ii) = cvec(ii+n1stop)
440 : end do
441 :
442 : ! *** ******************************************************************
443 : ! *** SECOND PART: Recomposition of the total wave scattered by
444 : ! *** the entire agglomerate/cluster by adding the
445 : ! *** waves scattered by each monomer taking into
446 : ! *** account the respective phase of the single waves.
447 : ! *** ******************************************************************
448 :
449 : ! *** ******************************************************************
450 : ! ----------------------------------------------------------------------
451 : ! 1) Calculate the amplitude functions |S1^j(th)| et |S2^j(th)|
452 : ! of one monomer of the agglomerate/cluster:
453 : ! ( see e.g. Bohren, Huffman (1983) p.112, Eq.(4.74) with the
454 : ! substitutions a_n -> d^1_1,n and b_n -> d^2_1,n
455 : ! or Rannou (1999) Eq.(1)-(6) )
456 : ! ----------------------------------------------------------------------
457 : ! *** ******************************************************************
458 :
459 0 : do iy=0,ntheta-1,1 ! loop over angles
460 0 : angle=iy*180._f/(ntheta-1)
461 0 : s1(iy)=0._f
462 0 : s2(iy)=0._f
463 0 : wom=cos(angle*3.1415926353_f/180._f)
464 :
465 0 : do ir=1,n1stop ! loop over Mie - indices
466 0 : xfact=2._f*(2._f*ir+1._f)/(ir*1._f*(ir*1._f+1._f))
467 0 : ajt=d1(ir)*fpi(ir,wom,fx_vars)+d2(ir)*tau(ir,wom,fx_vars)
468 0 : s1(iy)=s1(iy)+xfact*ajt
469 0 : ajt=d1(ir)*tau(ir,wom,fx_vars)+d2(ir)*fpi(ir,wom,fx_vars)
470 0 : s2(iy)=s2(iy)+xfact*ajt
471 : end do
472 :
473 0 : s11(iy)=abs(s1(iy))**2._f+abs(s2(iy))**2._f
474 0 : pol(iy)=abs(s1(iy))**2._f-abs(s2(iy))**2._f
475 0 : pol(iy)=pol(iy)/(abs(s1(iy))**2_f+abs(s2(iy))**2._f)
476 : ! *** S_11(theta) = 1/2 * ( |S_1|^2 + |S_2|^2 )
477 : ! *** above, s1(theta) = 2 * S_1(theta)
478 : ! *** =>S_11(theta) = 1/2 * ( |1/2*s1|^2 + |1/2*s2|^2 )
479 : ! = 1/8 * ( |s1|^2 + |s2|^2 )
480 0 : s11_n(iy)=.125_f*(abs(s1(iy))**2._f+abs(s2(iy))**2._f)
481 : end do
482 :
483 0 : s01=s1(0)
484 0 : s02=s2(0)
485 :
486 : ! *** Extinction cross section sigext( d^1_1,n , d^2_1,n ) ***
487 0 : sigext=0._f
488 0 : do ir=1,n1stop ! loop (sum) over Mie-indices
489 0 : sigext=sigext+(2._f*ir+1._f)*REAL(d1(ir)+d2(ir),f)
490 : end do
491 0 : sigext = fx_vars%nb * 2._f*pi/fx_vars%k**2._f * sigext ! Eq.(27)
492 :
493 : ! *** Alternatively (in a test, all values agreed with rel.acc. 1e-6),
494 : ! *** Extinction cross section sigext( S(0 deg) ) (optical theorem) ***
495 : ! *** (see e.g. Bohren, Huffman (1983), Eq. (4.76))
496 : ! *** S(0)=S_1(0)=S_2(0); sigma_ext = 4 pi / k^2 * Re(S(0))
497 : ! *** above, s1(theta) = 2 * S_1(theta) (factor 2 in 'xfact')
498 : ! sigext2 = nb * 4._f*pi/k**2._f * 0.5_f*REAL(s01)
499 : ! sigext3 = nb * 4._f*pi/k**2._f * 0.5_f*REAL(s02)
500 : ! *** ******************************************************************
501 :
502 : ! *** ******************************************************************
503 : ! ----------------------------------------------------------------------
504 : ! 2) Calculate the phase integral in Eq.(26) with P(r) already
505 : ! substituted ( compare Eq.(10) and (Botet 1995) ) :
506 : ! INT(0;infinity)[ sin(2XuZ) u^(d-2) f_co(u) du ]
507 : ! taking into account the different phases of the single
508 : ! scattered waves.
509 : ! ----------------------------------------------------------------------
510 : ! *** ******************************************************************
511 0 : do q=0,ntheta-1,1
512 0 : angle=q*180._f/(ntheta-1)
513 0 : if (angle .eq. 0._f) angle=0.001_f
514 0 : if (angle .eq. 180._f) angle=179.999_f
515 0 : fx_vars%zed=sin(angle*3.1415928353_f/180._f/2._f)
516 :
517 0 : bound=0._f
518 0 : interv=1
519 0 : errrel=1e-5_f
520 0 : p1=0._f
521 0 : dp1=0._f
522 :
523 : !======================================================================
524 : !--- Version using the QUADPACK - routine :
525 : !----------------------------------------------------------------------
526 0 : ifail = 0
527 0 : CALL dqagi(fp1,fx_vars,bound,interv,errabs,errrel,p1,dp1,neval,ifail,maxsub,lenw,nsubin,iwork,work)
528 0 : if(ifail.ne.0) then
529 0 : if (carma%f_do_print) write(carma%f_LUNOPRT, *) "fractal_meanfield_mod::ifail=",ifail," returned by dqagi()!"
530 0 : rc = RC_ERROR
531 0 : return
532 : endif
533 : !======================================================================
534 :
535 0 : p1=2._f*pi * (fx_vars%nb-1._f) / (fx_vars%coeff*fx_vars%zed*krg)*p1 + 1._f
536 0 : xint(q)=p1
537 : end do
538 :
539 : ! *** now, xint(theta) contains the square bracket terms in
540 : ! *** Botet (1997) Eq.(26) or Rannou (1999) Eq.(1)
541 :
542 : ! *** ******************************************************************
543 : ! ----------------------------------------------------------------------
544 : ! 3) Calculation of the phase function, calculation of the optical
545 : ! properties (asymmetrie factor g, scatt. cross section sigma_s)
546 : ! by angular integration: INT_0^180[ ... d_theta ]
547 : ! ----------------------------------------------------------------------
548 : ! *** ******************************************************************
549 :
550 : total=0._f
551 : totg=0._f
552 :
553 0 : do q=1,ntheta-2,2
554 0 : angle=(q-1)*180._f/(ntheta-1) ! angle in deg
555 0 : a0=fx_vars%nb*xint(q-1)*s11(q-1)*sin(angle*3.1415926353_f/180._f)
556 0 : c0=cos(angle*3.1415926353_f/180._f)
557 :
558 0 : angle=q*180._f/(ntheta-1)
559 0 : a1=fx_vars%nb*xint(q)*s11(q)*sin(angle*3.1415926353_f/180._f)
560 0 : c1=cos(angle*3.1415926353_f/180._f)
561 :
562 0 : angle=(q+1)*180._f/(ntheta-1)
563 0 : a2=fx_vars%nb*xint(q+1)*s11(q+1)*sin(angle*3.1415926353_f/180._f)
564 0 : c2=cos(angle*3.1415926353_f/180._f)
565 :
566 0 : total=total+2._f/6._f*3.1415926353_f/(ntheta-1)*(a0+4._f*a1+a2)
567 0 : totg=totg+2._f/6._f*3.1415926353_f/(ntheta-1)*(a0*c0+4._f*a1*c1+a2*c2)
568 : end do
569 0 : totg=totg/total
570 :
571 : ! *** ******************************************************************
572 : ! *** angular integration of I(theta) according to
573 : ! *** Botet (1997) Eq.(26) or Rannou (1999) Eq.(1)
574 : ! *** I(theta) = N 2pi/k^2 * S(theta) * [ phase integral ]
575 : ! *** with
576 : ! *** S(theta) = s11_n(i)
577 : ! *** [ phase i. ] = xint(i)
578 : ! *** Perfom integration using the following rule:
579 : ! *** Integral_0^pi[ I(theta) sin(theta) d_theta ]
580 : ! ***
581 : ! *** = Sum_q=1^ntheta-1{ Integral_th_(i-1)^th_i [
582 : ! ***
583 : ! *** 1/2(I(th_(i-1))+I(th_i)) * sin(th) d_th ] }
584 : ! ***
585 : ! *** = sin(delta_theta/2) * Sum_q=1^ntheta-1{
586 : ! ***
587 : ! *** ( I(th_(i-1)) + I(th_i) ) * sin(th_middle) }
588 : ! ***
589 : ! *** ******************************************************************
590 :
591 : !dthetad = 180._f / (ntheta-1) ! angular interval in deg
592 0 : dthetar = pi / (ntheta-1) ! angular interval in rad
593 0 : sumsca = 0._f
594 0 : do q=1,ntheta-1,1
595 0 : angler = (DBLE(q)-.5_f)*dthetar ! middle of interval in rad
596 0 : weight = SIN(angler) ! integration weight
597 0 : val = s11_n(q-1)*xint(q-1) + s11_n(q)*xint(q)
598 0 : sumsca = sumsca + val*weight
599 : end do
600 :
601 0 : sumsca = sin(.5_f*dthetar) * sumsca ! interval width factor
602 : ! *** Scattering cross section
603 0 : sigsca = 2._f * pi / fx_vars%k**2._f * DBLE(fx_vars%nb) * sumsca
604 : ! Warning! sigabs is well computed using this approximation
605 0 : sigabs=fx_vars%nb*(sigmae-sigmas)
606 : ! sigext=sigabs+sigsca is better than the mean-field value
607 : ! previously defined. This is used hereafter. (P.Rannou)
608 :
609 : ! *** Radius of equivalent mass sphere
610 0 : rems = fx_vars%a * fx_vars%nb**(1._f/3._f)
611 :
612 : ! *** reference area in definition of efficiencies is the geometrical
613 : ! *** cross section of equivalent mass sphere
614 0 : gems = pi * rems**2._f
615 :
616 : ! *** Extinction and scattering efficiencies:
617 0 : qsca = sigsca / gems
618 0 : qabs = sigabs / gems
619 0 : qext = qabs + qsca
620 :
621 0 : gfac = totg
622 :
623 : end subroutine fractal_meanfield
624 :
625 : !!
626 : !! Mie-scattering routine calling interface
627 : !!
628 : !! @author P. Rannou, R. Botet, Eric Wolf
629 : !! @version March 2013
630 0 : subroutine cmie(lambda,xn,xk,rad,an,bn,nstop)
631 :
632 : ! Arguments
633 : real(kind=f), intent(in) :: lambda !! wavelength (microns)
634 : real(kind=f), intent(in) :: xn !! real index of refraction
635 : real(kind=f), intent(in) :: xk !! imaginary index of refraction
636 : real(kind=f), intent(in) :: rad !! monomer radius (meters)
637 : complex(kind=f), intent(out) :: an(50) !! Mie wave coefficient an
638 : complex(kind=f), intent(out) :: bn(50) !! Mie wave coefficient bn
639 : integer, intent(out) :: nstop !! index of last mie-coefficent
640 :
641 : ! Local declarations
642 : integer, parameter :: nang = 451 ! number of angles
643 : complex(kind=f) :: refrel ! complex index of refraction
644 : real(kind=f) :: theta(10000)
645 : real(kind=f) :: x,dang
646 :
647 0 : refrel=cmplx(xn,xk,kind=f)
648 0 : x=2._f*3.14159265_f*rad/lambda ! size parameter of monomer
649 0 : dang=1.570796327_f/real(nang-1,kind=f)
650 :
651 0 : call intmie(x,refrel,nang,an,bn,nstop)
652 :
653 0 : return
654 : end subroutine cmie
655 :
656 : !!
657 : !! Mie scattering calculations
658 : !!
659 : !! @author P. Rannou, R. Botet, Eric Wolf
660 : !! @version March 2013
661 0 : SUBROUTINE intmie(x,refrel,nang,an,bn,nstop)
662 :
663 : ! Arguments
664 : real(kind=f), intent(in) :: x !! size parameter of monomer
665 : complex(kind=f), intent(in) :: refrel !! complex index of refraction
666 : integer, intent(in) :: nang !! number of angles
667 : complex(kind=f), intent(out) :: an(nf) !! Mie wave coefficient an
668 : complex(kind=f), intent(out) :: bn(nf) !! Mie wave coefficient an
669 : integer, intent(out) :: nstop !! index of last mie-coefficent
670 :
671 : ! Local declarations
672 : real(kind=f) :: amu(10000),pi(10000)
673 : real(kind=f) :: pi0(10000),pi1(10000)
674 : complex(kind=f) :: d(3000),y,xi,xi0,xi1
675 : complex(kind=f) :: s1(2000),s2(2000)
676 : real(kind=f) psi0,psi1,psi,dn,dx
677 : integer :: nmx,nn,n,j
678 : real(kind=f) :: rn, xstop, dang, ymod, chi0, chi1, apsi0, apsi1, fn, chi, apsi
679 :
680 0 : dx=x
681 0 : y=x*refrel
682 :
683 0 : xstop=x+4._f*x**.3333_f+2._f
684 0 : nstop=xstop
685 0 : ymod=abs(y)
686 0 : nmx=dmax1(xstop,ymod)+15
687 0 : dang=1.570796327_f/real(nang-1,kind=f)
688 :
689 : ! Initializations
690 0 : pi0(:) = 0._f
691 0 : pi1(:) = 0._f
692 0 : s1(:) = cmplx(0._f,0._f,kind=f)
693 0 : s2(:) = cmplx(0._f,0._f,kind=f)
694 0 : amu(:) = 0.0_f
695 0 : pi(:) = 0.0_f
696 :
697 0 : d(:) = cmplx(0._f,0._f,kind=f)
698 0 : nn=nmx-1
699 :
700 0 : do n=1,nn
701 0 : rn=nmx-n+1
702 0 : d(nmx-n)=(rn/y)-(1._f/(d(nmx-n+1)+rn/y))
703 : end do
704 :
705 0 : do j=1,nang
706 0 : pi0(j)=0._f ! Legendre functions
707 0 : pi1(j)=1._f
708 : end do
709 :
710 0 : nn=2*nang-1
711 :
712 0 : do j=1,nn
713 0 : s1(j)=cmplx(0._f,0._f,kind=f)
714 0 : s2(j)=cmplx(0._f,0._f,kind=f)
715 : end do
716 :
717 0 : psi0=cos(dx) ! Initialize Bessel functions
718 0 : psi1=sin(dx)
719 0 : chi0=-sin(x)
720 0 : chi1=cos(x)
721 :
722 0 : apsi0=psi0
723 0 : apsi1=psi1
724 :
725 0 : xi0=cmplx(apsi0,-chi0,kind=f)
726 0 : xi1=cmplx(apsi1,-chi1,kind=f)
727 :
728 0 : n=1
729 :
730 : ! ************* iterate over index n *************
731 0 : 200 dn=n
732 0 : rn=n
733 0 : fn=(2._f*rn+1._f)/(rn*(rn+1._f))
734 :
735 0 : psi=(2._f*dn-1._f)*psi1/dx-psi0 ! calculate Bessel functions
736 0 : chi=(2._f*rn-1._f)*chi1/x-chi0
737 0 : apsi=psi
738 0 : xi=cmplx(apsi,-chi,kind=f)
739 :
740 0 : an(n)=(d(n)/refrel+rn/x)*apsi-apsi1
741 0 : an(n)=an(n)/((d(n)/refrel+rn/x)*xi-xi1)
742 0 : bn(n)=(refrel*d(n)+rn/x)*apsi-apsi1
743 0 : bn(n)=bn(n)/((refrel*d(n)+rn/x)*xi-xi1)
744 :
745 0 : psi0=psi1
746 0 : psi1=psi
747 0 : apsi1=psi1
748 :
749 0 : chi0=chi1
750 0 : chi1=chi
751 0 : xi1=cmplx(apsi1,-chi1,kind=f)
752 :
753 0 : n=n+1
754 0 : rn=n
755 :
756 0 : do 999 j=1,nang
757 0 : pi1(j)=((2._f*rn-1._f)/(rn-1._f))*amu(j)*pi(j)
758 0 : pi1(j)=pi1(j)-rn*pi0(j)/(rn-1._f)
759 0 : 999 pi0(j)=pi(j)
760 :
761 0 : if (n-1-nstop) 200,300,300
762 : 300 continue
763 :
764 0 : return
765 : END SUBROUTINE intmie
766 :
767 : !!
768 : !!
769 : !! CALLS: FUNCTION dqag/dqdag/DADAPT_() Integration
770 : !! FUNCTION fpl() Integrand
771 : !!
772 : !! Integral in eq. 14, Botet et al. 1997
773 : !!
774 : !! @author P. Rannou, R. Botet, Eric Wolf
775 : !! @version March 2013
776 0 : function funa(carma,fx_vars,nu,n,p,rc)
777 : type(carma_type), intent(in) :: carma !! the carma object
778 : type(adgaquad_vars_type), intent(inout) :: fx_vars !! varaibles for functions being integrated
779 : integer, intent(in) :: n !! indices
780 : integer, intent(in) :: nu !! indices
781 : integer, intent(in) :: p !! indices
782 : integer, intent(inout) :: rc !! return code
783 : real(kind=f) :: funa !!
784 :
785 : ! Local declarations
786 : integer, parameter :: maxsub=1000
787 : real(kind=f) :: r,xa,xb,era,erl
788 : integer :: interv
789 : integer :: ifail
790 : integer, parameter :: lenw=4000 ! .ge. 4*maxsub
791 : integer :: iwork(maxsub),neval,nsubin ! nsubin=last
792 : real(kind=f) :: work(lenw)
793 : real(kind=f) :: bound, rres, rerr
794 :
795 : ! Set return code assuming success.
796 0 : rc = RC_OK
797 :
798 : ! Initializations
799 0 : funa=0._f
800 0 : fx_vars%u1=n
801 0 : fx_vars%u2=1
802 0 : fx_vars%u3=nu
803 0 : fx_vars%u4=1
804 0 : fx_vars%u5=p
805 0 : fx_vars%u6=0
806 0 : xa=-1._f
807 0 : xb=1._f
808 0 : bound=0._f
809 0 : interv=1
810 0 : era=0._f
811 0 : erl=1.e-4_f
812 0 : rres=0._f
813 0 : rerr=0._f
814 :
815 : !======================================================================
816 : !--- Version using the QUADPACK - routine :
817 : !----------------------------------------------------------------------
818 0 : ifail = 0
819 :
820 0 : call dqag(fpl,fx_vars,xa,xb,era,erl,3,rres,rerr,neval,ifail,maxsub,lenw,nsubin,iwork,work)
821 :
822 0 : if (ifail.ne.0) then
823 0 : if (carma%f_do_print) then
824 0 : write(carma%f_LUNOPRT, *) "funa::ifail=",ifail, &
825 0 : " returned by dqag() during call of funa(",nu,",",n,",",p,")"
826 : end if
827 0 : rc = RC_ERROR
828 0 : return
829 : endif
830 :
831 0 : rres=rres-2._f ! ceci est un artifice pour eviter que
832 : ! la routine se plante quand la fonction
833 : ! est paire (res=0.;err=1.d-3 impossible
834 : ! a atteindre!! j'ai fpl'=fpl+1....d'ou
835 : ! int(fpl)=int(fpl')-2. integr de -1 a 1!
836 :
837 0 : r = (2._f*p+1._f)/2._f
838 0 : funa = r * rres
839 :
840 0 : return
841 : END FUNCTION funa
842 :
843 : !!
844 : !! CALLS: FUNCTION calcpoly() Legendre-Functions
845 : !!
846 : !! Used in funa. Integrand of eq. 14, Botet et al. 1997
847 : !!
848 : !! @author P. Rannou, R. Botet, Eric Wolf
849 : !! @version March 2013
850 0 : FUNCTION fpl(x, fx_vars)
851 :
852 : ! Arguments
853 : real(kind=f),intent(in) :: x !!
854 : type(adgaquad_vars_type), intent(inout) :: fx_vars !! varaibles for functions being integrated
855 :
856 : ! Local declarations
857 : real(kind=f) :: fpl
858 : integer :: m,n,mu,nu,p,pmu
859 : real(kind=f) :: c1,c2,c3
860 :
861 0 : c1 = calcpoly(fx_vars%u1,fx_vars%u2,x,fx_vars)
862 0 : c2 = calcpoly(fx_vars%u3,fx_vars%u4,x,fx_vars)
863 0 : c3 = calcpoly(fx_vars%u5,fx_vars%u6,x,fx_vars)
864 :
865 0 : fpl=c1*c2*c3+1._f !this is a trick!
866 :
867 : return
868 : END FUNCTION fpl
869 :
870 : !!
871 : !! Calculate Legendre Polynomials, used in eq. 14 Botet et al. 1997
872 : !!
873 : !! @author P. Rannou, R. Botet, Eric Wolf
874 : !! @version March 2013
875 0 : FUNCTION calcpoly(l,m,x,fx_vars)
876 :
877 : ! Arguments
878 : integer, intent(in) :: l !! indices
879 : integer, intent(in) :: m !! indices
880 : real(kind=f), intent(in) :: x !! return result
881 : type(adgaquad_vars_type), intent(in) :: fx_vars !! variables for functions being integrated
882 :
883 : ! Local declarations
884 : real(kind=f) :: calcpoly
885 : integer ::lbl
886 : real(kind=f) :: pll, pmm, somx2, pmmp1
887 : integer :: i, ll
888 : real(kind=f) :: fact1
889 : integer :: mstar
890 :
891 0 : mstar=m
892 :
893 0 : lbl=0
894 0 : calcpoly=0._f
895 :
896 0 : if (mstar.lt.0)then
897 0 : mstar=-m
898 0 : lbl=1
899 : endif
900 :
901 0 : if (mstar.gt.l) then
902 0 : pll=0._f
903 0 : calcpoly=0._f
904 : return ! si m>l, Pl,m=0 !
905 : endif
906 :
907 0 : pmm=1._f
908 :
909 0 : if(mstar.gt.0) then
910 0 : somx2=sqrt((1._f-x)*(1._f+x))
911 0 : fact1=1._f
912 0 : do i=1,mstar
913 0 : pmm=+pmm*fact1*somx2 !cghmt - en + !!
914 0 : fact1=fact1+2._f
915 : end do
916 : endif
917 :
918 0 : if(l.eq.mstar) then
919 : calcpoly=pmm
920 : else
921 0 : pmmp1=x*(2*mstar+1)*pmm
922 :
923 0 : if(l.eq.mstar+1) then
924 : calcpoly=pmmp1
925 : else
926 0 : do ll=mstar+2,l
927 0 : pll=(x*(2*ll-1)*pmmp1-(ll+mstar-1)*pmm)/(ll-mstar)
928 0 : pmm=pmmp1
929 0 : pmmp1=pll
930 : end do
931 : calcpoly=pll
932 : endif
933 : endif
934 :
935 0 : if (lbl.eq.1) then
936 0 : calcpoly=(-1)**mstar*(fx_vars%fact(l-mstar)/fx_vars%fact(l+mstar))*calcpoly
937 0 : mstar=-m !restitution du parametre m!!!!!
938 : endif
939 :
940 : return
941 : END FUNCTION calcpoly
942 :
943 : !!
944 : !! replaces funb(nu,n,p) in original code,
945 : !! saving n*n re-calculations of funa(nu,n,p).
946 : !!
947 : !! Calculates eq. 15, Botet et al. 1997
948 : !!
949 : !! @author P. Rannou, R. Botet, Eric Wolf
950 : !! @version March 2013
951 0 : FUNCTION funb_n(nu,n,p,funca)
952 :
953 : ! Arguments
954 : integer, intent(in) :: nu !! indices
955 : integer, intent(in) :: n !! indices
956 : integer, intent(in) :: p !! indices
957 : real(kind=f), intent(in) :: funca(nmi,nmi,0:n2m) !! return result
958 :
959 : ! Local Declarations
960 : real(kind=f) :: funb_n
961 : integer :: i, l, j
962 : real(kind=f) :: var
963 :
964 0 : funb_n = 0._f
965 0 : i = int((p*1._f-1._f-abs(n*1._f-nu*1._f))*1._f/2._f)
966 : !print*,nu,n,p,i
967 :
968 0 : do l=0,i
969 0 : j = p-2*l-1
970 :
971 : ! omit j = -1 (when nu=n and p=l=i=0)
972 0 : IF (j .GE. 0) THEN
973 :
974 0 : var = funca(nu,n,j) ! in main, a(nu,n,p) was stored in
975 : ! funca(nu,n;p)
976 0 : funb_n = funb_n + var
977 : ENDIF
978 :
979 : end do
980 0 : funb_n = (2._f*p+1._f) * funb_n
981 : return
982 : END FUNCTION funb_n
983 :
984 : !! Replaces funs(pp,k) in original code
985 : !!
986 : !! CALLS:
987 : !! FUNCTION dqagi/dq2agi/DADAPT_() Integration
988 : !! FUNCTION xfreal_n() Integrand
989 : !! FUNCTION xfimag_n() Integrand
990 : !!
991 : !! Calculates eq. 16 , Botet et al. 1997
992 : !!
993 : !! @author P. Rannou, R. Botet, Eric Wolf
994 : !! @version Mar 2013
995 0 : function funs_n(carma,fx_vars,pp,rc)
996 :
997 : ! Arguments
998 : type(carma_type), intent(in) :: carma !! the carma object
999 : type(adgaquad_vars_type), intent(inout) :: fx_vars !! varaibles for functions being integrated
1000 : integer, intent(in) :: pp !! indices
1001 : integer, intent(inout) :: rc !! return code
1002 :
1003 : ! Local Declarations
1004 : integer, parameter :: maxsub=50000
1005 : complex(kind=f) :: rcomplex,funs_n
1006 : real(kind=f) :: rres,ires,rerr,ierr,afun
1007 : real(kind=f) :: xa,xb
1008 : integer :: ifail
1009 : integer, parameter :: lenw=200000 ! .ge. 4*maxsub
1010 : integer :: iwork(maxsub),neval,nsubin ! nsubin=last
1011 : real(kind=f) :: work(lenw)
1012 : real(kind=f) :: rg, bound, errabs, errrel
1013 : integer :: interv
1014 :
1015 0 : rc = RC_OK
1016 :
1017 0 : rg=fx_vars%alpha*fx_vars%nb**(1._f/fx_vars%df)*fx_vars%a
1018 :
1019 0 : afun=(2._f*3.1415926_f)/(fx_vars%k**3._f)
1020 :
1021 :
1022 0 : fx_vars%pbes=pp
1023 0 : fx_vars%kbes=fx_vars%k
1024 :
1025 0 : bound=0._f
1026 0 : interv=1
1027 :
1028 0 : errabs=0._f
1029 0 : errrel=1.e-3_f
1030 :
1031 0 : rres=0._f
1032 : !trres=0._f
1033 0 : rerr=0._f
1034 : !trerr=0._f
1035 0 : xa=0._f
1036 0 : xb=5._f*fx_vars%k*rg
1037 :
1038 : !======================================================================
1039 : !--- Version using the QUADPACK - routine :
1040 : !----------------------------------------------------------------------
1041 0 : ifail = 0
1042 0 : CALL dqagi(xfreal_n,fx_vars,bound,interv,errabs,errrel,rres,rerr,neval,ifail,maxsub,lenw,nsubin,iwork,work)
1043 0 : if (ifail.ne.0) then
1044 0 : if (carma%f_do_print) then
1045 0 : write(carma%f_LUNOPRT, *) "funs_n::ifail=",ifail, &
1046 0 : " returned by dqag() during call of funs(",pp, &
1047 0 : ") while integrating xfreal_n()"
1048 : end if
1049 0 : rc = RC_ERROR
1050 0 : return
1051 : endif
1052 :
1053 0 : bound=0._f
1054 0 : interv=1
1055 :
1056 0 : ires=0._f
1057 0 : ierr=0._f
1058 0 : xa=0._f
1059 0 : xb=5._f*fx_vars%k*rg
1060 :
1061 : !======================================================================
1062 : !--- Version using the QUADPACK - routine :
1063 : !----------------------------------------------------------------------
1064 0 : ifail = 0
1065 0 : CALL dqagi(xfimag_n,fx_vars,bound,interv,errabs,errrel,ires,ierr,neval,ifail,maxsub,lenw,nsubin,iwork,work)
1066 0 : if(ifail.ne.0) then
1067 0 : if (carma%f_do_print) then
1068 0 : write(carma%f_LUNOPRT, *) "funs_n::ifail=",ifail, &
1069 0 : " returned by dqagi() during call of funs(",pp, &
1070 0 : ") while integrating xfimag_n()"
1071 : end if
1072 0 : rc = RC_ERROR
1073 0 : return
1074 : endif
1075 :
1076 0 : rcomplex = cmplx(1._f,0._f,kind=f)*rres + cmplx(0._f,1._f,kind=f)*ires
1077 :
1078 0 : funs_n = afun * rcomplex
1079 :
1080 : continue
1081 0 : return
1082 : END FUNCTION funs_n
1083 :
1084 : !!
1085 : !! replaces xfreel(xx) in original code
1086 : !! CALLS: FUNCTION BESSELJY() Spherical Bessel functions
1087 : !! FUNCTION phi() Probability distrib.
1088 : !!
1089 : !! @author P. Rannou, R. Botet, Eric Wolf
1090 : !! @version Mar 2013
1091 0 : FUNCTION xfreal_n(xx, fx_vars)
1092 :
1093 : ! Arguments
1094 : real(kind=f), intent(in) :: xx
1095 : type(adgaquad_vars_type), intent(inout) :: fx_vars !! varaibles for functions being integrated
1096 :
1097 : ! Local Declarations
1098 : complex(kind=f) :: z,xj(0:nf),xjp(0:nf),xy(0:nf),xyp(0:nf)
1099 : complex(kind=f) :: jsol,ysol,hsol,hpsol
1100 : real(kind=f) :: x,r,xfreal_n
1101 : integer :: ifail,p,pc
1102 : real(kind=f) :: rg
1103 :
1104 0 : ifail = 0
1105 0 : rg = fx_vars%alpha*fx_vars%nb**(1._f/fx_vars%df)*fx_vars%a
1106 0 : x = xx
1107 0 : if (x.GT.3000._f) x=3000._f
1108 0 : z = x*cmplx(1._f,0._f,kind=f)
1109 :
1110 0 : pc = fx_vars%pbes
1111 0 : if( fx_vars%pbes .eq. 0 ) pc = fx_vars%pbes + 1
1112 :
1113 0 : CALL BESSELJY(z,pc,xj,xjp,xy,xyp,ifail)
1114 :
1115 0 : r=x/fx_vars%kbes
1116 :
1117 0 : xfreal_n = real( z*z*xj(fx_vars%pbes)*xj(fx_vars%pbes) * phi(r,fx_vars) )
1118 :
1119 : return
1120 : END FUNCTION xfreal_n
1121 :
1122 : !!
1123 : !! replaces xfima(xx) in original code
1124 : !! CALLS: FUNCTION BESSELJY() Spherical Bessel functions
1125 : !!
1126 : !! @author P. Rannou, R. Botet, Eric Wolf
1127 : !! @version Mar 2013
1128 0 : FUNCTION xfimag_n(xx, fx_vars)
1129 :
1130 : ! Arguments
1131 : real(kind=f), intent(in) :: xx
1132 : type(adgaquad_vars_type), intent(inout) :: fx_vars !! variables for functions being integrated
1133 :
1134 : ! Local Declarations
1135 : complex(kind=f) :: z,xj(0:nf),xjp(0:nf),xy(0:nf),xyp(0:nf)
1136 : real(kind=f) :: x,r,xfimag_n
1137 : integer :: ifail,p,pc
1138 : real(kind=f) :: rg
1139 :
1140 0 : rg = fx_vars%alpha*fx_vars%nb**(1._f/fx_vars%df)*fx_vars%a
1141 0 : x = xx
1142 0 : if (x.gt.3000._f) x=3000._f
1143 0 : ifail = 0
1144 0 : z = x*cmplx(1._f,0._f,kind=f)
1145 :
1146 0 : pc = fx_vars%pbes
1147 0 : if( fx_vars%pbes .eq. 0 ) pc = fx_vars%pbes + 1
1148 :
1149 0 : CALL BESSELJY(z,pc,xj,xjp,xy,xyp,ifail)
1150 :
1151 0 : r=x/fx_vars%kbes
1152 :
1153 0 : xfimag_n = real( z*z*xj(fx_vars%pbes)*xy(fx_vars%pbes) * phi(r,fx_vars) )
1154 :
1155 : return
1156 : END FUNCTION xfimag_n
1157 :
1158 :
1159 : !! Spherical Bessel functions j_n(z) and y_n(z) of complex
1160 : !! argument to desired accuracy,
1161 : !! and their derivatives, up to a maximal order n=LMAX.
1162 : !! j_n(z) = SQRT(pi/2 / z) * J_(n + 1/2)(z)
1163 : !! y_n(z) = SQRT(pi/2 / z) * Y_(n + 1/2)(z)
1164 : !! Adapted from:
1165 : !! I.J.Thompson, A.R.Barnett
1166 : !! "Modified Bessel Funkctions I_v(z) and K_v(z)
1167 : !! of Real Order and Complex Argument, to Selected
1168 : !! Accuracy"
1169 : !! COMP.PHYS.COMMUN. 47 (1987) 245-57
1170 : !! (Source code printed on page 249)
1171 : !! ******************************************************************
1172 : !! INPUTS:
1173 : !! X argument z, dble cmplx
1174 : !! z in the upper half plane, Im(z) > -3
1175 : !! LMAX largest desired order of Bessel functions, int
1176 : !! j_n,y_n,j_n',y_n' are calculated for n=0 to n=LMAX
1177 : !! Dimension of arrays xj,xjp,xy,xyp at least (0:LMAX)
1178 : !! XJ(M) Spher. Bessel function j_m(z), dble cmplx
1179 : !! XJP(M) Derivative of Spher. Bessel function d/dz [ j_m(z) ],
1180 : !! dble cmplx
1181 : !! XY(M) Spher. Bessel function y_m(z), dble cmplx
1182 : !! XYP(M) Derivative of Spher. Bessel function d/dz [ y_m(z) ],
1183 : !! dble cmplx
1184 : !! IFAIL error flag, int
1185 : !! = 0 if all results are satisfactory
1186 : !! = -1 for arguments out of range
1187 : !! = > 0 for results ok up to and including the
1188 : !! function of order LMAX-IFAIL
1189 : !!
1190 : !! @author P. Rannou, R. Botet, Eric Wolf
1191 : !! @version Mar 2013
1192 0 : SUBROUTINE BESSELJY (X, LMAX, XJ, XJP, XY, XYP, IFAIL)
1193 :
1194 : ! Arguments
1195 : complex(kind=f), intent(in) :: X
1196 : integer, intent(in) :: LMAX
1197 : complex(kind=f), intent(out) :: XJ(0:LMAX)
1198 : complex(kind=f), intent(out) :: XJP(0:LMAX)
1199 : complex(kind=f), intent(out) :: XY(0:LMAX)
1200 : complex(kind=f), intent(out) :: XYP(0:LMAX)
1201 : integer, intent(out) :: IFAIL
1202 :
1203 : ! Local Declarations
1204 : INTEGER, PARAMETER :: LIMIT = 20000
1205 : REAL(kind=f),parameter :: ZERO = 0._f
1206 : REAL(kind=f),parameter :: ONE = 1._f
1207 : REAL(kind=f),parameter :: ACCUR = 1e-12_f
1208 : REAL(kind=f),parameter :: TM30 = 1e-30_f
1209 : COMPLEX(kind=f), parameter :: CI = (0._f, 1._f)
1210 : complex(kind=f) :: XI, W, PL, B, D, FF, DEL, C, XJ0, XH1, XH1P, XTEMP
1211 : integer :: L
1212 :
1213 0 : IF (ABS(X).LT.ACCUR .OR. AIMAG(X) .LT. -3.d0) THEN
1214 0 : IFAIL=-1
1215 0 : GOTO 5
1216 : END IF
1217 :
1218 : ! *** Lentz - Algorithmus (?) :
1219 0 : XI = ONE/X
1220 0 : W = XI + XI
1221 0 : PL = LMAX * XI
1222 0 : FF = PL + XI
1223 0 : B = FF + FF + XI
1224 0 : D = ZERO
1225 0 : C = FF
1226 0 : DO 1 L=1,LIMIT
1227 0 : D = B - D
1228 0 : C = B - ONE/C
1229 0 : IF(ABS(D).LT. TM30) D = TM30
1230 0 : IF(ABS(C).LT. TM30) C = TM30
1231 0 : D = ONE / D
1232 0 : DEL = D * C
1233 0 : FF = FF * DEL
1234 0 : B = B + W
1235 0 : 1 IF(ABS(DEL-ONE).LT.ACCUR) GOTO 2
1236 0 : IFAIL = -2
1237 0 : GOTO 5
1238 :
1239 0 : 2 XJ(LMAX) = TM30
1240 0 : XJP(LMAX) = FF * XJ(LMAX)
1241 :
1242 : ! *** Abwaertsrekursion
1243 0 : DO 3 L = LMAX-1,0,-1
1244 0 : XJ(L) = PL * XJ(L+1) + XJP(L+1)
1245 0 : XJP(L) = PL * XJ(L) - XJ(L+1)
1246 0 : 3 PL = PL - XI
1247 :
1248 : ! *** Calculate the l=0 Besselfunktionen
1249 0 : XJ0 = XI * SIN(X)
1250 0 : XY(0) = - XI * COS(X)
1251 0 : XH1 = EXP(CI * X) * XI * (-CI)
1252 0 : XH1P = XH1 * (CI - XI)
1253 0 : B = XH1P
1254 :
1255 : ! *** Rescale XJ, XJP, converting to spherical Bessels
1256 : ! *** Recur XH1,XH1P as sperical Bessels
1257 0 : W = ONE / XJ(0)
1258 0 : PL = XI
1259 0 : DO 4 L = 0,LMAX
1260 0 : XJ(L) = XJ0 * (W*XJ(L))
1261 0 : XJP(L) = XJ0 * (W*XJP(L)) - XI * XJ(L)
1262 0 : IF (L.EQ.0) GOTO 4
1263 0 : XTEMP = XH1
1264 0 : XH1 = (PL-XI) * XTEMP - XH1P
1265 0 : PL = PL + XI
1266 0 : XH1P = - PL * XH1 + XTEMP
1267 0 : XY(L) = CI * (XJ(L) - XH1) ! y_n = i * ( j_n - h^1_n )
1268 0 : XYP(L) = CI * (XJP(L) - XH1P) ! und dito fuer Ableitungen
1269 0 : 4 CONTINUE
1270 0 : XYP(0) = CI * (XJP(0) - B)
1271 0 : RETURN
1272 :
1273 0 : 5 WRITE(*,10) IFAIL
1274 : 10 FORMAT( 'ERROR in SUBR BESSELJY() : IFAIL = ', I4)
1275 0 : RETURN
1276 : END SUBROUTINE BESSELJY
1277 :
1278 : !!
1279 : !! Angular function pi_l( x=cos(theta) )
1280 : !! e.g. Bohren,Huffman (1983)
1281 : !! pp.94 ff Eq.(4.46)-(4.49)
1282 : !! p.112
1283 : !! CALLS: FUNCTION calcpoly() Legendre-Functions
1284 : !!
1285 : !! @author P. Rannou, R. Botet, Eric Wolf
1286 : !! @version Mar 2013
1287 0 : FUNCTION fpi(l,x,fx_vars)
1288 :
1289 : ! Arguments
1290 : integer, intent(in) :: l
1291 : real(kind=f), intent(in) :: x
1292 : type(adgaquad_vars_type), intent(inout) :: fx_vars !! varaibles for functions being integrated
1293 :
1294 : ! Local declarations
1295 : real(kind=f) :: fpi
1296 : real(kind=f) :: y
1297 : real(kind=f) :: flag
1298 :
1299 0 : y=x
1300 0 : if (x.eq.1._f) y=1._f-1.e-6_f
1301 : ! alternatively, one could use Bohren,Huffman
1302 : ! p.112: pi_n(1)=tau_n(1)= 1/2 * n * (n+1) !!!
1303 0 : flag=calcpoly(l,1,y,fx_vars)
1304 0 : fpi=(1._f-y**2._f)**(-0.5_f)*flag
1305 : return
1306 : END FUNCTION fpi
1307 :
1308 : !!
1309 : !! Angular function tau_l( x=cos(theta) )
1310 : !! e.g. Bohren,Huffman (1983)
1311 : !! pp.94 ff Eq.(4.46)-(4.49)
1312 : !! p.112
1313 : !! CALLS: FUNCTION calcpoly() Legendre-Functions
1314 : !!
1315 : !! @author P. Rannou, R. Botet, Eric Wolf
1316 : !! @version March 2013
1317 0 : FUNCTION tau(l,x,fx_vars)
1318 :
1319 : ! Arguments
1320 : integer, intent(in) :: l
1321 : real(kind=f), intent(in) :: x
1322 : type(adgaquad_vars_type), intent(inout) :: fx_vars !! varaibles for functions being integrated
1323 :
1324 : ! Local Declarations
1325 : real(kind=f) :: fp
1326 : real(kind=f) :: tau
1327 : real(kind=f) :: flag
1328 : real(kind=f) :: y
1329 :
1330 0 : y=x
1331 0 : if (x.eq.1._f) y=1._f-1.e-6_f
1332 : ! alternatively, one could use Bohren,Huffman
1333 : ! p.112: pi_n(1)=tau_n(1)= 1/2 * n * (n+1) !!!
1334 0 : flag=calcpoly(l,0,y,fx_vars)
1335 0 : fp=fpi(l,y,fx_vars)
1336 0 : tau=-y*fp+l*(l*1._f+1._f)*flag
1337 : return
1338 : END FUNCTION tau
1339 :
1340 : !!
1341 : !!
1342 : !!
1343 : !! @author P. Rannou, R. Botet, Eric Wolf
1344 : !! @version March 2013
1345 0 : function fp1(u, fx_vars)
1346 :
1347 : ! Arguments
1348 : real(kind=f), intent(in) :: u !!
1349 : type(adgaquad_vars_type), intent(inout) :: fx_vars !! varaibles for functions being integrated
1350 : real(kind=f) :: fp1 !! returns
1351 :
1352 : ! Local Declarations
1353 : real(kind=f) :: krg,s1,s2,s3,rg
1354 :
1355 0 : rg=fx_vars%alpha*fx_vars%a*fx_vars%nb**(1._f/fx_vars%df)
1356 0 : krg=fx_vars%k*rg
1357 0 : s1=sin(2._f*krg*fx_vars%zed*u)
1358 0 : s2=u**(fx_vars%df-2._f)
1359 0 : s3=fco(u, fx_vars)
1360 0 : fp1=s1*s2*s3
1361 :
1362 : return
1363 : END FUNCTION fp1
1364 :
1365 : !!
1366 : !! CALLS: FUNCTION dqagi/dqdagi/DADAPT_() Integration
1367 : !! FUNCTION fdval() Integrand
1368 : !!
1369 : !! @author P. Rannou, R. Botet, Eric Wolf
1370 : !! @version March 2013
1371 0 : FUNCTION anorm(carma, fx_vars, rc)
1372 :
1373 : ! arguments
1374 : type(carma_type), intent(in) :: carma !! the carma object
1375 : type(adgaquad_vars_type), intent(inout) :: fx_vars !! varaibles for functions being integrated
1376 : integer, intent(inout) :: rc !! return code
1377 :
1378 : ! Local Declarations
1379 : real(kind=f) :: anorm
1380 : integer :: interv
1381 : integer, parameter :: maxsub=50000
1382 : integer :: ifail
1383 : integer, parameter :: lenw=200000 ! .ge. 4*maxsub
1384 : integer :: iwork(maxsub),neval,nsubin ! nsubin=last
1385 : real(kind=f) :: work(lenw)
1386 : real(kind=f) :: bound,errrel,errabs,b,db,c
1387 :
1388 0 : rc = RC_OK
1389 :
1390 0 : bound=0._f
1391 0 : interv=1
1392 0 : errrel=1.e-3_f
1393 0 : errabs=0._f
1394 0 : b=0._f
1395 0 : db=0._f
1396 :
1397 : !======================================================================
1398 : !--- Version using the QUADPACK - routine :
1399 : !----------------------------------------------------------------------
1400 0 : ifail = 0
1401 0 : CALL dqagi(fdval,fx_vars,bound,interv,errabs,errrel,b,db,neval,ifail,maxsub,lenw,nsubin,iwork,work)
1402 0 : if(ifail.ne.0) then
1403 0 : if (carma%f_do_print) write(carma%f_LUNOPRT, *) "anorm::ifail=",ifail," returned by dqagi() during call of anorm"
1404 0 : rc = RC_ERROR
1405 0 : return
1406 : endif
1407 :
1408 0 : c=0.5_f
1409 0 : anorm=b*4._f*3.1415926_f
1410 0 : return
1411 : END FUNCTION anorm
1412 :
1413 : !!
1414 : !! Probability distribution of monomer location within cluster
1415 : !! CALLS: FUNCTION fdval()
1416 : !!
1417 : !! @author P. Rannou, R. Botet, Eric Wolf
1418 : !! @version March 2013
1419 0 : FUNCTION phi(x,fx_vars)
1420 :
1421 : ! Arguments
1422 : real(kind=f), intent(in) :: x
1423 : type(adgaquad_vars_type), intent(inout) :: fx_vars !! varaibles for functions being integrated
1424 :
1425 : ! Local Declarations
1426 : real(kind=f) :: fval,pref,phi, rg, z
1427 :
1428 0 : rg=fx_vars%alpha*fx_vars%nb**(1._f/fx_vars%df)*fx_vars%a
1429 0 : z=x/rg
1430 0 : pref=(x/rg)**(fx_vars%df-3._f)/(fx_vars%coeff*rg**3._f)
1431 0 : fval=z**(1._f-fx_vars%df)*fdval(z, fx_vars)
1432 0 : phi=pref*fval
1433 : continue
1434 : return
1435 : END FUNCTION phi
1436 :
1437 : !!
1438 : !! Probability distribution of monomer location within cluster
1439 : !! CALLS: FUNCTION fdval()
1440 : !!
1441 : !! @author P. Rannou, R. Botet, Eric Wolf
1442 : !! @version March 2013
1443 0 : FUNCTION fco(z, fx_vars)
1444 :
1445 : ! Arguments
1446 : real(kind=f), intent(in) :: z
1447 : type(adgaquad_vars_type), intent(inout) :: fx_vars !! varaibles for functions being integrated
1448 :
1449 : ! Local Declarations
1450 : real(kind=f) :: fco
1451 : real(kind=f) :: fval
1452 :
1453 0 : fval=z**(1._f-fx_vars%df)*fdval(z, fx_vars)
1454 0 : fco=fval
1455 : continue
1456 : return
1457 : END FUNCTION fco
1458 :
1459 : !!
1460 : !! @author P. Rannou, R. Botet, Eric Wolf
1461 : !! @version March 2013
1462 0 : FUNCTION fdval(x, fx_vars)
1463 :
1464 : type(adgaquad_vars_type), intent(inout) :: fx_vars !! varaibles for functions being integrated
1465 :
1466 : ! Arguments
1467 : real(kind=f), intent(in) :: x
1468 :
1469 : ! Local Declarations
1470 : real(kind=f) :: fdval
1471 :
1472 0 : fdval=x**(fx_vars%df-1._f)*exp(-x**fx_vars%df/2._f)
1473 : return
1474 : END FUNCTION fdval
1475 :
1476 : end module fractal_meanfield_mod
|
