subroutine miev0 ( xx, crefin, perfct, mimcut, anyang, 1
$ numang, xmu, nmom, ipolzn, momdim, prnt,
$ qext, qsca, gqsc, pmom, sforw, sback, s1,
$ s2, tforw, tback )
c
c computes mie scattering and extinction efficiencies; asymmetry
c factor; forward- and backscatter amplitude; scattering
c amplitudes for incident polarization parallel and perpendicular
c to the plane of scattering, as functions of scattering angle;
c coefficients in the legendre polynomial expansions of either the
c unpolarized phase function or the polarized phase matrix;
c and some quantities needed in polarized radiative transfer.
c
c calls : biga, ckinmi, small1, small2, testmi, miprnt,
c lpcoef, errmsg
c
c i n t e r n a l v a r i a b l e s
c -----------------------------------
c
c an,bn mie coefficients little-a-sub-n, little-b-sub-n
c ( ref. 1, eq. 16 )
c anm1,bnm1 mie coefficients little-a-sub-(n-1),
c little-b-sub-(n-1); used in -gqsc- sum
c anp coeffs. in s+ expansion ( ref. 2, p. 1507 )
c bnp coeffs. in s- expansion ( ref. 2, p. 1507 )
c anpm coeffs. in s+ expansion ( ref. 2, p. 1507 )
c when mu is replaced by - mu
c bnpm coeffs. in s- expansion ( ref. 2, p. 1507 )
c when mu is replaced by - mu
c calcmo(k) true, calculate moments for k-th phase quantity
c (derived from -ipolzn-; used only in 'lpcoef')
c cbiga(n) bessel function ratio capital-a-sub-n (ref. 2, eq. 2)
c ( complex version )
c cior complex index of refraction with negative
c imaginary part (van de hulst convention)
c cioriv 1 / cior
c coeff ( 2n + 1 ) / ( n ( n + 1 ) )
c fn floating point version of index in loop performing
c mie series summation
c lita,litb(n) mie coefficients -an-, -bn-, saved in arrays for
c use in calculating legendre moments *pmom*
c maxtrm max. possible no. of terms in mie series
c mm + 1 and - 1, alternately.
c mim magnitude of imaginary refractive index
c mre real part of refractive index
c maxang max. possible value of input variable -numang-
c nangd2 (numang+1)/2 ( no. of angles in 0-90 deg; anyang=f )
c noabs true, sphere non-absorbing (determined by -mimcut-)
c np1dn ( n + 1 ) / n
c npquan highest-numbered phase quantity for which moments are
c to be calculated (the largest digit in -ipolzn-
c if ipolzn .ne. 0)
c ntrm no. of terms in mie series
c pass1 true on first entry, false thereafter; for self-test
c pin(j) angular function little-pi-sub-n ( ref. 2, eq. 3 )
c at j-th angle
c pinm1(j) little-pi-sub-(n-1) ( see -pin- ) at j-th angle
c psinm1 ricatti-bessel function psi-sub-(n-1), argument -xx-
c psin ricatti-bessel function psi-sub-n of argument -xx-
c ( ref. 1, p. 11 ff. )
c rbiga(n) bessel function ratio capital-a-sub-n (ref. 2, eq. 2)
c ( real version, for when imag refrac index = 0 )
c rioriv 1 / mre
c rn 1 / n
c rtmp (real) temporary variable
c sp(j) s+ for j-th angle ( ref. 2, p. 1507 )
c sm(j) s- for j-th angle ( ref. 2, p. 1507 )
c sps(j) s+ for (numang+1-j)-th angle ( anyang=false )
c sms(j) s- for (numang+1-j)-th angle ( anyang=false )
c taun angular function little-tau-sub-n ( ref. 2, eq. 4 )
c at j-th angle
c tcoef n ( n+1 ) ( 2n+1 ) (for summing tforw,tback series)
c twonp1 2n + 1
c yesang true if scattering amplitudes are to be calculated
c zetnm1 ricatti-bessel function zeta-sub-(n-1) of argument
c -xx- ( ref. 2, eq. 17 )
c zetn ricatti-bessel function zeta-sub-n of argument -xx-
c
c ----------------------------------------------------------------------
c -------- i / o specifications for subroutine miev0 -----------------
c ----------------------------------------------------------------------
implicit none
logical anyang, perfct, prnt(*)
integer ipolzn, momdim, numang, nmom
real*8 gqsc, mimcut, pmom( 0:momdim, * ), qext, qsca,
$ xmu(*), xx
complex*16 crefin, sforw, sback, s1(*), s2(*), tforw(*),
$ tback(*)
integer maxang,mxang2,maxtrm
real*8 onethr
c ----------------------------------------------------------------------
c
parameter ( maxang = 501, mxang2 = maxang/2 + 1 )
c
c ** note -- maxtrm = 10100 is neces-
c ** sary to do some of the test probs,
c ** but 1100 is sufficient for most
c ** conceivable applications
parameter ( maxtrm = 1100 )
parameter ( onethr = 1./3. )
c
logical anysav, calcmo(4), noabs, ok, pass1, persav, yesang
integer npquan
integer i,j,n,nmosav,iposav,numsav,ntrm,nangd2
real*8 mim, mimsav, mre, mm, np1dn
real*8 rioriv,xmusav,xxsav,sq,fn,rn,twonp1,tcoef, coeff
real*8 xinv,psinm1,chinm1,psin,chin,rtmp,taun
real*8 rbiga( maxtrm ), pin( maxang ), pinm1( maxang )
complex*16 an, bn, anm1, bnm1, anp, bnp, anpm, bnpm, cresav,
$ cior, cioriv, ctmp, zet, zetnm1, zetn
complex*16 cbiga( maxtrm ), lita( maxtrm ), litb( maxtrm ),
$ sp( maxang ), sm( maxang ), sps( mxang2 ), sms( mxang2 )
equivalence ( cbiga, rbiga )
save pass1
sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
data pass1 / .true. /
c
c
if ( pass1 ) then
c ** save certain user input values
xxsav = xx
cresav = crefin
mimsav = mimcut
persav = perfct
anysav = anyang
nmosav = nmom
iposav = ipolzn
numsav = numang
xmusav = xmu( 1 )
c ** reset input values for test case
xx = 10.0
crefin = ( 1.5, - 0.1 )
perfct = .false.
mimcut = 0.0
anyang = .true.
numang = 1
xmu( 1 )= - 0.7660444
nmom = 1
ipolzn = - 1
c
end if
c ** check input and calculate
c ** certain variables from input
c
10 call ckinmi
( numang, maxang, xx, perfct, crefin, momdim,
$ nmom, ipolzn, anyang, xmu, calcmo, npquan )
c
if ( perfct .and. xx .le. 0.1 ) then
c ** use totally-reflecting
c ** small-particle limit
c
call small1 ( xx, numang, xmu, qext, qsca, gqsc, sforw,
$ sback, s1, s2, tforw, tback, lita, litb )
ntrm = 2
go to 200
end if
c
if ( .not.perfct ) then
c
cior = crefin
if ( dimag( cior ) .gt. 0.0 ) cior = dconjg( cior )
mre = dble( cior )
mim = - dimag( cior )
noabs = mim .le. mimcut
cioriv = 1.0 / cior
rioriv = 1.0 / mre
c
if ( xx * dmax1( 1.d0, cdabs(cior) ) .le. 0.d1 ) then
c
c ** use general-refractive-index
c ** small-particle limit
c ** ( ref. 2, p. 1508 )
c
call small2 ( xx, cior, .not.noabs, numang, xmu, qext,
$ qsca, gqsc, sforw, sback, s1, s2, tforw,
$ tback, lita, litb )
ntrm = 2
go to 200
end if
c
end if
c
nangd2 = ( numang + 1 ) / 2
yesang = numang .gt. 0
c ** estimate number of terms in mie series
c ** ( ref. 2, p. 1508 )
if ( xx.le.8.0 ) then
ntrm = xx + 4. * xx**onethr + 1.
else if ( xx.lt.4200. ) then
ntrm = xx + 4.05 * xx**onethr + 2.
else
ntrm = xx + 4. * xx**onethr + 2.
end if
if ( ntrm+1 .gt. maxtrm )
$ call errmsg( 'miev0--parameter maxtrm too small', .true. )
c
c ** calculate logarithmic derivatives of
c ** j-bessel-fcn., big-a-sub-(1 to ntrm)
if ( .not.perfct )
$ call biga( cior, xx, ntrm, noabs, yesang, rbiga, cbiga )
c
c ** initialize ricatti-bessel functions
c ** (psi,chi,zeta)-sub-(0,1) for upward
c ** recurrence ( ref. 1, eq. 19 )
xinv = 1.0 / xx
psinm1 = dsin( xx )
chinm1 = dcos( xx )
psin = psinm1 * xinv - chinm1
chin = chinm1 * xinv + psinm1
zetnm1 = dcmplx( psinm1, chinm1 )
zetn = dcmplx( psin, chin )
c ** initialize previous coeffi-
c ** cients for -gqsc- series
anm1 = ( 0.0, 0.0 )
bnm1 = ( 0.0, 0.0 )
c ** initialize angular function little-pi
c ** and sums for s+, s- ( ref. 2, p. 1507 )
if ( anyang ) then
do 60 j = 1, numang
pinm1( j ) = 0.0
pin( j ) = 1.0
sp ( j ) = ( 0.0, 0.0 )
sm ( j ) = ( 0.0, 0.0 )
60 continue
else
do 70 j = 1, nangd2
pinm1( j ) = 0.0
pin( j ) = 1.0
sp ( j ) = ( 0.0, 0.0 )
sm ( j ) = ( 0.0, 0.0 )
sps( j ) = ( 0.0, 0.0 )
sms( j ) = ( 0.0, 0.0 )
70 continue
end if
c ** initialize mie sums for efficiencies, etc.
qsca = 0.0
gqsc = 0.0
sforw = ( 0., 0. )
sback = ( 0., 0. )
tforw( 1 ) = ( 0., 0. )
tback( 1 ) = ( 0., 0. )
c
c
c --------- loop to sum mie series -----------------------------------
c
mm = + 1.0
do 100 n = 1, ntrm
c ** compute various numerical coefficients
fn = n
rn = 1.0 / fn
np1dn = 1.0 + rn
twonp1 = 2 * n + 1
coeff = twonp1 / ( fn * ( n + 1 ) )
tcoef = twonp1 * ( fn * ( n + 1 ) )
c
c ** calculate mie series coefficients
if ( perfct ) then
c ** totally-reflecting case
c
an = ( ( fn*xinv ) * psin - psinm1 ) /
$ ( ( fn*xinv ) * zetn - zetnm1 )
bn = psin / zetn
c
else if ( noabs ) then
c ** no-absorption case
c
an = ( ( rioriv*rbiga(n) + ( fn*xinv ) ) * psin - psinm1 )
$ / ( ( rioriv*rbiga(n) + ( fn*xinv ) ) * zetn - zetnm1 )
bn = ( ( mre * rbiga(n) + ( fn*xinv ) ) * psin - psinm1 )
$ / ( ( mre * rbiga(n) + ( fn*xinv ) ) * zetn - zetnm1 )
else
c ** absorptive case
c
an = ( ( cioriv * cbiga(n) + ( fn*xinv ) ) * psin - psinm1 )
$ /( ( cioriv * cbiga(n) + ( fn*xinv ) ) * zetn - zetnm1 )
bn = ( ( cior * cbiga(n) + ( fn*xinv ) ) * psin - psinm1 )
$ /( ( cior * cbiga(n) + ( fn*xinv ) ) * zetn - zetnm1 )
qsca = qsca + twonp1 * ( sq( an ) + sq( bn ) )
c
end if
c ** save mie coefficients for *pmom* calculation
lita( n ) = an
litb( n ) = bn
c ** increment mie sums for non-angle-
c ** dependent quantities
c
sforw = sforw + twonp1 * ( an + bn )
tforw( 1 ) = tforw( 1 ) + tcoef * ( an - bn )
sback = sback + ( mm * twonp1 ) * ( an - bn )
tback( 1 ) = tback( 1 ) + ( mm * tcoef ) * ( an + bn )
gqsc = gqsc + ( fn - rn ) * dble( anm1 * dconjg( an )
$ + bnm1 * dconjg( bn ) )
$ + coeff * dble( an * dconjg( bn ) )
c
if ( yesang ) then
c ** put mie coefficients in form
c ** needed for computing s+, s-
c ** ( ref. 2, p. 1507 )
anp = coeff * ( an + bn )
bnp = coeff * ( an - bn )
c ** increment mie sums for s+, s-
c ** while upward recursing
c ** angular functions little pi
c ** and little tau
if ( anyang ) then
c ** arbitrary angles
c
c ** vectorizable loop
do 80 j = 1, numang
rtmp = ( xmu( j ) * pin( j ) ) - pinm1( j )
taun = fn * rtmp - pinm1( j )
sp( j ) = sp( j ) + anp * ( pin( j ) + taun )
sm( j ) = sm( j ) + bnp * ( pin( j ) - taun )
pinm1( j ) = pin( j )
pin( j ) = ( xmu( j ) * pin( j ) ) + np1dn * rtmp
80 continue
c
else
c ** angles symmetric about 90 degrees
anpm = mm * anp
bnpm = mm * bnp
c ** vectorizable loop
do 90 j = 1, nangd2
rtmp = ( xmu( j ) * pin( j ) ) - pinm1( j )
taun = fn * rtmp - pinm1( j )
sp ( j ) = sp ( j ) + anp * ( pin( j ) + taun )
sms( j ) = sms( j ) + bnpm * ( pin( j ) + taun )
sm ( j ) = sm ( j ) + bnp * ( pin( j ) - taun )
sps( j ) = sps( j ) + anpm * ( pin( j ) - taun )
pinm1( j ) = pin( j )
pin( j ) = ( xmu( j ) * pin( j ) ) + np1dn * rtmp
90 continue
c
end if
end if
c ** update relevant quantities for next
c ** pass through loop
mm = - mm
anm1 = an
bnm1 = bn
c ** upward recurrence for ricatti-bessel
c ** functions ( ref. 1, eq. 17 )
c
zet = ( twonp1 * xinv ) * zetn - zetnm1
zetnm1 = zetn
zetn = zet
psinm1 = psin
psin = dble( zetn )
100 continue
c
c ---------- end loop to sum mie series --------------------------------
c
c
qext = 2. / xx**2 * dble( sforw )
if ( perfct .or. noabs ) then
qsca = qext
else
qsca = 2. / xx**2 * qsca
end if
c
gqsc = 4. / xx**2 * gqsc
sforw = 0.5 * sforw
sback = 0.5 * sback
tforw( 2 ) = 0.5 * ( sforw + 0.25 * tforw( 1 ) )
tforw( 1 ) = 0.5 * ( sforw - 0.25 * tforw( 1 ) )
tback( 2 ) = 0.5 * ( sback + 0.25 * tback( 1 ) )
tback( 1 ) = 0.5 * ( - sback + 0.25 * tback( 1 ) )
c
if ( yesang ) then
c ** recover scattering amplitudes
c ** from s+, s- ( ref. 1, eq. 11 )
if ( anyang ) then
c ** vectorizable loop
do 110 j = 1, numang
s1( j ) = 0.5 * ( sp( j ) + sm( j ) )
s2( j ) = 0.5 * ( sp( j ) - sm( j ) )
110 continue
c
else
c ** vectorizable loop
do 120 j = 1, nangd2
s1( j ) = 0.5 * ( sp( j ) + sm( j ) )
s2( j ) = 0.5 * ( sp( j ) - sm( j ) )
120 continue
c ** vectorizable loop
do 130 j = 1, nangd2
s1( numang+1 - j ) = 0.5 * ( sps( j ) + sms( j ) )
s2( numang+1 - j ) = 0.5 * ( sps( j ) - sms( j ) )
130 continue
end if
c
end if
c ** calculate legendre moments
200 if ( nmom.gt.0 )
$ call lpcoef ( ntrm, nmom, ipolzn, momdim, calcmo, npquan,
$ lita, litb, pmom )
c
if ( dimag(crefin) .gt. 0.0 ) then
c ** take complex conjugates
c ** of scattering amplitudes
sforw = dconjg( sforw )
sback = dconjg( sback )
do 210 i = 1, 2
tforw( i ) = dconjg( tforw(i) )
tback( i ) = dconjg( tback(i) )
210 continue
c
do 220 j = 1, numang
s1( j ) = dconjg( s1(j) )
s2( j ) = dconjg( s2(j) )
220 continue
c
end if
c
if ( pass1 ) then
c ** compare test case results with
c ** correct answers and abort if bad
c
call testmi ( qext, qsca, gqsc, sforw, sback, s1, s2,
$ tforw, tback, pmom, momdim, ok )
if ( .not. ok ) then
prnt(1) = .false.
prnt(2) = .false.
call miprnt( prnt, xx, perfct, crefin, numang, xmu, qext,
$ qsca, gqsc, nmom, ipolzn, momdim, calcmo,
$ pmom, sforw, sback, tforw, tback, s1, s2 )
call errmsg( 'miev0 -- self-test failed', .true. )
end if
c ** restore user input values
xx = xxsav
crefin = cresav
mimcut = mimsav
perfct = persav
anyang = anysav
nmom = nmosav
ipolzn = iposav
numang = numsav
xmu(1) = xmusav
pass1 = .false.
go to 10
c
end if
c
if ( prnt(1) .or. prnt(2) )
$ call miprnt( prnt, xx, perfct, crefin, numang, xmu, qext,
$ qsca, gqsc, nmom, ipolzn, momdim, calcmo,
$ pmom, sforw, sback, tforw, tback, s1, s2 )
c
return
c
end
subroutine ckinmi( numang, maxang, xx, perfct, crefin, momdim, 1
$ nmom, ipolzn, anyang, xmu, calcmo, npquan )
c
c check for bad input to 'miev0' and calculate -calcmo,npquan-
c
implicit none
logical perfct, anyang, calcmo(*)
integer numang, maxang, momdim, nmom, ipolzn, npquan
real*8 xx, xmu(*)
integer i,l,j,ip
complex*16 crefin
c
character*4 string
logical inperr
c
inperr = .false.
c
if ( numang.gt.maxang ) then
call errmsg( 'miev0--parameter maxang too small', .true. )
inperr = .true.
end if
if ( numang.lt.0 ) call wrtbad( 'numang', inperr )
if ( xx.lt.0. ) call wrtbad( 'xx', inperr )
if ( .not.perfct .and. dble(crefin).le.0. )then
print *,'crefin=',crefin
call wrtbad( 'crefin', inperr )
endif
if ( momdim.lt.1 ) call wrtbad( 'momdim', inperr )
c
if ( nmom.ne.0 ) then
if ( nmom.lt.0 .or. nmom.gt.momdim ) call wrtbad('nmom',inperr)
if ( iabs(ipolzn).gt.4444 ) call wrtbad( 'ipolzn', inperr )
npquan = 0
do 5 l = 1, 4
calcmo( l ) = .false.
5 continue
if ( ipolzn.ne.0 ) then
c ** parse out -ipolzn- into its digits
c ** to find which phase quantities are
c ** to have their moments calculated
c
write( string, '(i4)' ) iabs(ipolzn)
do 10 j = 1, 4
ip = ichar( string(j:j) ) - ichar( '0' )
if ( ip.ge.1 .and. ip.le.4 ) calcmo( ip ) = .true.
if ( ip.eq.0 .or. (ip.ge.5 .and. ip.le.9) )
$ call wrtbad( 'ipolzn', inperr )
npquan = max0( npquan, ip )
10 continue
end if
end if
c
if ( anyang ) then
c ** allow for slight imperfections in
c ** computation of cosine
do 20 i = 1, numang
if ( xmu(i).lt.-1.00001 .or. xmu(i).gt.1.00001 )
$ call wrtbad( 'xmu', inperr )
20 continue
else
do 22 i = 1, ( numang + 1 ) / 2
if ( xmu(i).lt.-0.00001 .or. xmu(i).gt.1.00001 )
$ call wrtbad( 'xmu', inperr )
22 continue
end if
c
if ( inperr )
$ call errmsg( 'miev0--input error(s). aborting...', .true. )
c
if ( xx.gt.20000.0 .or. dble(crefin).gt.10.0 .or.
$ dabs( dimag(crefin) ).gt.10.0 ) call errmsg(
$ 'miev0--xx or crefin outside tested range', .false. )
c
return
end
subroutine lpcoef ( ntrm, nmom, ipolzn, momdim, calcmo, npquan, 1
$ a, b, pmom )
c
c calculate legendre polynomial expansion coefficients (also
c called moments) for phase quantities ( ref. 5 formulation )
c
c input: ntrm number terms in mie series
c nmom, ipolzn, momdim 'miev0' arguments
c calcmo flags calculated from -ipolzn-
c npquan defined in 'miev0'
c a, b mie series coefficients
c
c output: pmom legendre moments ('miev0' argument)
c
c *** notes ***
c
c (1) eqs. 2-5 are in error in dave, appl. opt. 9,
c 1888 (1970). eq. 2 refers to m1, not m2; eq. 3 refers to
c m2, not m1. in eqs. 4 and 5, the subscripts on the second
c term in square brackets should be interchanged.
c
c (2) the general-case logic in this subroutine works correctly
c in the two-term mie series case, but subroutine 'lpco2t'
c is called instead, for speed.
c
c (3) subroutine 'lpco1t', to do the one-term case, is never
c called within the context of 'miev0', but is included for
c complete generality.
c
c (4) some improvement in speed is obtainable by combining the
c 310- and 410-loops, if moments for both the third and fourth
c phase quantities are desired, because the third phase quantity
c is the real part of a complex series, while the fourth phase
c quantity is the imaginary part of that very same series. but
c most users are not interested in the fourth phase quantity,
c which is related to circular polarization, so the present
c scheme is usually more efficient.
c
implicit none
logical calcmo(*)
integer ipolzn, momdim, nmom, ntrm, npquan
real*8 pmom( 0:momdim, * )
complex*16 a(*), b(*)
c
c ** specification of local variables
c
c am(m) numerical coefficients a-sub-m-super-l
c in dave, eqs. 1-15, as simplified in ref. 5.
c
c bi(i) numerical coefficients b-sub-i-super-l
c in dave, eqs. 1-15, as simplified in ref. 5.
c
c bidel(i) 1/2 bi(i) times factor capital-del in dave
c
c cm,dm() arrays c and d in dave, eqs. 16-17 (mueller form),
c calculated using recurrence derived in ref. 5
c
c cs,ds() arrays c and d in ref. 4, eqs. a5-a6 (sekera form),
c calculated using recurrence derived in ref. 5
c
c c,d() either -cm,dm- or -cs,ds-, depending on -ipolzn-
c
c evenl true for even-numbered moments; false otherwise
c
c idel 1 + little-del in dave
c
c maxtrm max. no. of terms in mie series
c
c maxmom max. no. of non-zero moments
c
c nummom number of non-zero moments
c
c recip(k) 1 / k
c
integer maxtrm,maxmom,mxmom2,maxrcp
parameter ( maxtrm = 1102, maxmom = 2*maxtrm, mxmom2 = maxmom/2,
$ maxrcp = 4*maxtrm + 2 )
real*8 am( 0:maxtrm ), bi( 0:mxmom2 ), bidel( 0:mxmom2 ),
$ recip( maxrcp )
complex*16 cm( maxtrm ), dm( maxtrm ), cs( maxtrm ), ds( maxtrm ),
$ c( maxtrm ), d( maxtrm )
integer k,j,l,nummom,ld2,idel,m,i,mmax,imax
real*8 sum
equivalence ( c, cm ), ( d, dm )
logical pass1, evenl
save pass1, recip
data pass1 / .true. /
c
c
if ( pass1 ) then
c
do 1 k = 1, maxrcp
recip( k ) = 1.0 / k
1 continue
pass1 = .false.
c
end if
c
do 5 j = 1, max0( 1, npquan )
do 5 l = 0, nmom
pmom( l, j ) = 0.0
5 continue
c
if ( ntrm.eq.1 ) then
call lpco1t ( nmom, ipolzn, momdim, calcmo, a, b, pmom )
return
else if ( ntrm.eq.2 ) then
call lpco2t ( nmom, ipolzn, momdim, calcmo, a, b, pmom )
return
end if
c
if ( ntrm+2 .gt. maxtrm )
$ call errmsg( 'lpcoef--parameter maxtrm too small', .true. )
c
c ** calculate mueller c, d arrays
cm( ntrm+2 ) = ( 0., 0. )
dm( ntrm+2 ) = ( 0., 0. )
cm( ntrm+1 ) = ( 1. - recip( ntrm+1 ) ) * b( ntrm )
dm( ntrm+1 ) = ( 1. - recip( ntrm+1 ) ) * a( ntrm )
cm( ntrm ) = ( recip(ntrm) + recip(ntrm+1) ) * a( ntrm )
$ + ( 1. - recip(ntrm) ) * b( ntrm-1 )
dm( ntrm ) = ( recip(ntrm) + recip(ntrm+1) ) * b( ntrm )
$ + ( 1. - recip(ntrm) ) * a( ntrm-1 )
c
do 10 k = ntrm-1, 2, -1
cm( k ) = cm( k+2 ) - ( 1. + recip(k+1) ) * b( k+1 )
$ + ( recip(k) + recip(k+1) ) * a( k )
$ + ( 1. - recip(k) ) * b( k-1 )
dm( k ) = dm( k+2 ) - ( 1. + recip(k+1) ) * a( k+1 )
$ + ( recip(k) + recip(k+1) ) * b( k )
$ + ( 1. - recip(k) ) * a( k-1 )
10 continue
cm( 1 ) = cm( 3 ) + 1.5 * ( a( 1 ) - b( 2 ) )
dm( 1 ) = dm( 3 ) + 1.5 * ( b( 1 ) - a( 2 ) )
c
if ( ipolzn.ge.0 ) then
c
do 20 k = 1, ntrm + 2
c( k ) = ( 2*k - 1 ) * cm( k )
d( k ) = ( 2*k - 1 ) * dm( k )
20 continue
c
else
c ** compute sekera c and d arrays
cs( ntrm+2 ) = ( 0., 0. )
ds( ntrm+2 ) = ( 0., 0. )
cs( ntrm+1 ) = ( 0., 0. )
ds( ntrm+1 ) = ( 0., 0. )
c
do 30 k = ntrm, 1, -1
cs( k ) = cs( k+2 ) + ( 2*k + 1 ) * ( cm( k+1 ) - b( k ) )
ds( k ) = ds( k+2 ) + ( 2*k + 1 ) * ( dm( k+1 ) - a( k ) )
30 continue
c
do 40 k = 1, ntrm + 2
c( k ) = ( 2*k - 1 ) * cs( k )
d( k ) = ( 2*k - 1 ) * ds( k )
40 continue
c
end if
c
c
if( ipolzn.lt.0 ) nummom = min0( nmom, 2*ntrm - 2 )
if( ipolzn.ge.0 ) nummom = min0( nmom, 2*ntrm )
if ( nummom .gt. maxmom )
$ call errmsg( 'lpcoef--parameter maxtrm too small', .true. )
c
c ** loop over moments
do 500 l = 0, nummom
ld2 = l / 2
evenl = mod( l,2 ) .eq. 0
c ** calculate numerical coefficients
c ** a-sub-m and b-sub-i in dave
c ** double-sums for moments
if( l.eq.0 ) then
c
idel = 1
do 60 m = 0, ntrm
am( m ) = 2.0 * recip( 2*m + 1 )
60 continue
bi( 0 ) = 1.0
c
else if( evenl ) then
c
idel = 1
do 70 m = ld2, ntrm
am( m ) = ( 1. + recip( 2*m-l+1 ) ) * am( m )
70 continue
do 75 i = 0, ld2-1
bi( i ) = ( 1. - recip( l-2*i ) ) * bi( i )
75 continue
bi( ld2 ) = ( 2. - recip( l ) ) * bi( ld2-1 )
c
else
c
idel = 2
do 80 m = ld2, ntrm
am( m ) = ( 1. - recip( 2*m+l+2 ) ) * am( m )
80 continue
do 85 i = 0, ld2
bi( i ) = ( 1. - recip( l+2*i+1 ) ) * bi( i )
85 continue
c
end if
c ** establish upper limits for sums
c ** and incorporate factor capital-
c ** del into b-sub-i
mmax = ntrm - idel
if( ipolzn.ge.0 ) mmax = mmax + 1
imax = min0( ld2, mmax - ld2 )
if( imax.lt.0 ) go to 600
do 90 i = 0, imax
bidel( i ) = bi( i )
90 continue
if( evenl ) bidel( 0 ) = 0.5 * bidel( 0 )
c
c ** perform double sums just for
c ** phase quantities desired by user
if( ipolzn.eq.0 ) then
c
do 110 i = 0, imax
c ** vectorizable loop (cray)
sum = 0.0
do 100 m = ld2, mmax - i
sum = sum + am( m ) *
$ ( dble( c(m-i+1) * dconjg( c(m+i+idel) ) )
$ + dble( d(m-i+1) * dconjg( d(m+i+idel) ) ) )
100 continue
pmom( l,1 ) = pmom( l,1 ) + bidel( i ) * sum
110 continue
pmom( l,1 ) = 0.5 * pmom( l,1 )
go to 500
c
end if
c
if ( calcmo(1) ) then
do 160 i = 0, imax
c ** vectorizable loop (cray)
sum = 0.0
do 150 m = ld2, mmax - i
sum = sum + am( m ) *
$ dble( c(m-i+1) * dconjg( c(m+i+idel) ) )
150 continue
pmom( l,1 ) = pmom( l,1 ) + bidel( i ) * sum
160 continue
end if
c
c
if ( calcmo(2) ) then
do 210 i = 0, imax
c ** vectorizable loop (cray)
sum = 0.0
do 200 m = ld2, mmax - i
sum = sum + am( m ) *
$ dble( d(m-i+1) * dconjg( d(m+i+idel) ) )
200 continue
pmom( l,2 ) = pmom( l,2 ) + bidel( i ) * sum
210 continue
end if
c
c
if ( calcmo(3) ) then
do 310 i = 0, imax
c ** vectorizable loop (cray)
sum = 0.0
do 300 m = ld2, mmax - i
sum = sum + am( m ) *
$ ( dble( c(m-i+1) * dconjg( d(m+i+idel) ) )
$ + dble( c(m+i+idel) * dconjg( d(m-i+1) ) ) )
300 continue
pmom( l,3 ) = pmom( l,3 ) + bidel( i ) * sum
310 continue
pmom( l,3 ) = 0.5 * pmom( l,3 )
end if
c
c
if ( calcmo(4) ) then
do 410 i = 0, imax
c ** vectorizable loop (cray)
sum = 0.0
do 400 m = ld2, mmax - i
sum = sum + am( m ) *
$ ( dimag( c(m-i+1) * dconjg( d(m+i+idel) ) )
$ + dimag( c(m+i+idel) * dconjg( d(m-i+1) ) ))
400 continue
pmom( l,4 ) = pmom( l,4 ) + bidel( i ) * sum
410 continue
pmom( l,4 ) = - 0.5 * pmom( l,4 )
end if
c
500 continue
c
c
600 return
end
subroutine lpco1t ( nmom, ipolzn, momdim, calcmo, a, b, pmom ) 1
c
c calculate legendre polynomial expansion coefficients (also
c called moments) for phase quantities in special case where
c no. terms in mie series = 1
c
c input: nmom, ipolzn, momdim 'miev0' arguments
c calcmo flags calculated from -ipolzn-
c a(1), b(1) mie series coefficients
c
c output: pmom legendre moments
c
implicit none
logical calcmo(*)
integer ipolzn, momdim, nmom,nummom,l
real*8 pmom( 0:momdim, * ),sq,a1sq,b1sq
complex*16 a(*), b(*), ctmp, a1b1c
sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
c
c
a1sq = sq( a(1) )
b1sq = sq( b(1) )
a1b1c = a(1) * dconjg( b(1) )
c
if( ipolzn.lt.0 ) then
c
if( calcmo(1) ) pmom( 0,1 ) = 2.25 * b1sq
if( calcmo(2) ) pmom( 0,2 ) = 2.25 * a1sq
if( calcmo(3) ) pmom( 0,3 ) = 2.25 * dble( a1b1c )
if( calcmo(4) ) pmom( 0,4 ) = 2.25 *dimag( a1b1c )
c
else
c
nummom = min0( nmom, 2 )
c ** loop over moments
do 100 l = 0, nummom
c
if( ipolzn.eq.0 ) then
if( l.eq.0 ) pmom( l,1 ) = 1.5 * ( a1sq + b1sq )
if( l.eq.1 ) pmom( l,1 ) = 1.5 * dble( a1b1c )
if( l.eq.2 ) pmom( l,1 ) = 0.15 * ( a1sq + b1sq )
go to 100
end if
c
if( calcmo(1) ) then
if( l.eq.0 ) pmom( l,1 ) = 2.25 * ( a1sq + b1sq / 3. )
if( l.eq.1 ) pmom( l,1 ) = 1.5 * dble( a1b1c )
if( l.eq.2 ) pmom( l,1 ) = 0.3 * b1sq
end if
c
if( calcmo(2) ) then
if( l.eq.0 ) pmom( l,2 ) = 2.25 * ( b1sq + a1sq / 3. )
if( l.eq.1 ) pmom( l,2 ) = 1.5 * dble( a1b1c )
if( l.eq.2 ) pmom( l,2 ) = 0.3 * a1sq
end if
c
if( calcmo(3) ) then
if( l.eq.0 ) pmom( l,3 ) = 3.0 * dble( a1b1c )
if( l.eq.1 ) pmom( l,3 ) = 0.75 * ( a1sq + b1sq )
if( l.eq.2 ) pmom( l,3 ) = 0.3 * dble( a1b1c )
end if
c
if( calcmo(4) ) then
if( l.eq.0 ) pmom( l,4 ) = - 1.5 * dimag( a1b1c )
if( l.eq.1 ) pmom( l,4 ) = 0.0
if( l.eq.2 ) pmom( l,4 ) = 0.3 * dimag( a1b1c )
end if
c
100 continue
c
end if
c
return
end
subroutine lpco2t ( nmom, ipolzn, momdim, calcmo, a, b, pmom ) 1
c
c calculate legendre polynomial expansion coefficients (also
c called moments) for phase quantities in special case where
c no. terms in mie series = 2
c
c input: nmom, ipolzn, momdim 'miev0' arguments
c calcmo flags calculated from -ipolzn-
c a(1-2), b(1-2) mie series coefficients
c
c output: pmom legendre moments
c
implicit none
logical calcmo(*)
integer ipolzn, momdim, nmom,l,nummom
real*8 pmom( 0:momdim, * ),sq,a1sq,b1sq,pm1,pm2,a2sq,b2sq
complex*16 a(*), b(*)
complex*16 a2c, b2c, ctmp, ca, cac, cat, cb, cbc, cbt, cg, ch
sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
c
c
ca = 3. * a(1) - 5. * b(2)
cat= 3. * b(1) - 5. * a(2)
cac = dconjg( ca )
a2sq = sq( a(2) )
b2sq = sq( b(2) )
a2c = dconjg( a(2) )
b2c = dconjg( b(2) )
c
if( ipolzn.lt.0 ) then
c ** loop over sekera moments
nummom = min0( nmom, 2 )
do 50 l = 0, nummom
c
if( calcmo(1) ) then
if( l.eq.0 ) pmom( l,1 ) = 0.25 * ( sq(cat) +
$ (100./3.) * b2sq )
if( l.eq.1 ) pmom( l,1 ) = (5./3.) * dble( cat * b2c )
if( l.eq.2 ) pmom( l,1 ) = (10./3.) * b2sq
end if
c
if( calcmo(2) ) then
if( l.eq.0 ) pmom( l,2 ) = 0.25 * ( sq(ca) +
$ (100./3.) * a2sq )
if( l.eq.1 ) pmom( l,2 ) = (5./3.) * dble( ca * a2c )
if( l.eq.2 ) pmom( l,2 ) = (10./3.) * a2sq
end if
c
if( calcmo(3) ) then
if( l.eq.0 ) pmom( l,3 ) = 0.25 * dble( cat*cac +
$ (100./3.)*b(2)*a2c )
if( l.eq.1 ) pmom( l,3 ) = 5./6. * dble( b(2)*cac +
$ cat*a2c )
if( l.eq.2 ) pmom( l,3 ) = 10./3. * dble( b(2) * a2c )
end if
c
if( calcmo(4) ) then
if( l.eq.0 ) pmom( l,4 ) = -0.25 * dimag( cat*cac +
$ (100./3.)*b(2)*a2c )
if( l.eq.1 ) pmom( l,4 ) = -5./6. * dimag( b(2)*cac +
$ cat*a2c )
if( l.eq.2 ) pmom( l,4 ) = -10./3. * dimag( b(2) * a2c )
end if
c
50 continue
c
else
c
cb = 3. * b(1) + 5. * a(2)
cbt= 3. * a(1) + 5. * b(2)
cbc = dconjg( cb )
cg = ( cbc*cbt + 10.*( cac*a(2) + b2c*cat) ) / 3.
ch = 2.*( cbc*a(2) + b2c*cbt )
c
c ** loop over mueller moments
nummom = min0( nmom, 4 )
do 100 l = 0, nummom
c
if( ipolzn.eq.0 .or. calcmo(1) ) then
if( l.eq.0 ) pm1 = 0.25 * sq(ca) + sq(cb) / 12.
$ + (5./3.) * dble(ca*b2c) + 5.*b2sq
if( l.eq.1 ) pm1 = dble( cb * ( cac/6. + b2c ) )
if( l.eq.2 ) pm1 = sq(cb)/30. + (20./7.) * b2sq
$ + (2./3.) * dble( ca * b2c )
if( l.eq.3 ) pm1 = (2./7.) * dble( cb * b2c )
if( l.eq.4 ) pm1 = (40./63.) * b2sq
if ( calcmo(1) ) pmom( l,1 ) = pm1
end if
c
if( ipolzn.eq.0 .or. calcmo(2) ) then
if( l.eq.0 ) pm2 = 0.25*sq(cat) + sq(cbt) / 12.
$ + (5./3.) * dble(cat*a2c) + 5.*a2sq
if( l.eq.1 ) pm2 = dble( cbt * ( dconjg(cat)/6. + a2c) )
if( l.eq.2 ) pm2 = sq(cbt)/30. + (20./7.) * a2sq
$ + (2./3.) * dble( cat * a2c )
if( l.eq.3 ) pm2 = (2./7.) * dble( cbt * a2c )
if( l.eq.4 ) pm2 = (40./63.) * a2sq
if ( calcmo(2) ) pmom( l,2 ) = pm2
end if
c
if( ipolzn.eq.0 ) then
pmom( l,1 ) = 0.5 * ( pm1 + pm2 )
go to 100
end if
c
if( calcmo(3) ) then
if( l.eq.0 ) pmom( l,3 ) = 0.25 * dble( cac*cat + cg +
$ 20.*b2c*a(2) )
if( l.eq.1 ) pmom( l,3 ) = dble( cac*cbt + cbc*cat +
$ 3.*ch ) / 12.
if( l.eq.2 ) pmom( l,3 ) = 0.1 * dble( cg + (200./7.) *
$ b2c * a(2) )
if( l.eq.3 ) pmom( l,3 ) = dble( ch ) / 14.
if( l.eq.4 ) pmom( l,3 ) = 40./63. * dble( b2c * a(2) )
end if
c
if( calcmo(4) ) then
if( l.eq.0 ) pmom( l,4 ) = 0.25 * dimag( cac*cat + cg +
$ 20.*b2c*a(2) )
if( l.eq.1 ) pmom( l,4 ) = dimag( cac*cbt + cbc*cat +
$ 3.*ch ) / 12.
if( l.eq.2 ) pmom( l,4 ) = 0.1 * dimag( cg + (200./7.) *
$ b2c * a(2) )
if( l.eq.3 ) pmom( l,4 ) = dimag( ch ) / 14.
if( l.eq.4 ) pmom( l,4 ) = 40./63. * dimag( b2c * a(2) )
end if
c
100 continue
c
end if
c
return
end
subroutine biga( cior, xx, ntrm, noabs, yesang, rbiga, cbiga ) 1
c
c calculate logarithmic derivatives of j-bessel-function
c
c input : cior, xx, ntrm, noabs, yesang (defined in 'miev0')
c
c output : rbiga or cbiga (defined in 'miev0')
c
c internal variables :
c
c confra value of lentz continued fraction for -cbiga(ntrm)-,
c used to initialize downward recurrence.
c down = true, use down-recurrence. false, do not.
c f1,f2,f3 arithmetic statement functions used in determining
c whether to use up- or down-recurrence
c ( ref. 2, eqs. 6-8 )
c mre real refractive index
c mim imaginary refractive index
c rezinv 1 / ( mre * xx ); temporary variable for recurrence
c zinv 1 / ( cior * xx ); temporary variable for recurrence
c
implicit none
logical down, noabs, yesang
integer ntrm,n
real*8 mre, mim, rbiga(*), xx, rezinv, rtmp, f1,f2,f3
complex*16 cior, ctmp, confra, cbiga(*), zinv
f1( mre ) = - 8.0 + mre**2 * ( 26.22 + mre * ( - 0.4474
$ + mre**3 * ( 0.00204 - 0.000175 * mre ) ) )
f2( mre ) = 3.9 + mre * ( - 10.8 + 13.78 * mre )
f3( mre ) = - 15.04 + mre * ( 8.42 + 16.35 * mre )
c
c ** decide whether 'biga' can be
c ** calculated by up-recurrence
mre = dble( cior )
mim = dabs( dimag( cior ) )
if ( mre.lt.1.0 .or. mre.gt.10.0 .or. mim.gt.10.0 ) then
down = .true.
else if ( yesang ) then
down = .true.
if ( mim*xx .lt. f2( mre ) ) down = .false.
else
down = .true.
if ( mim*xx .lt. f1( mre ) ) down = .false.
end if
c
zinv = 1.0 / ( cior * xx )
rezinv = 1.0 / ( mre * xx )
c
if ( down ) then
c ** compute initial high-order 'biga' using
c ** lentz method ( ref. 1, pp. 17-20 )
c
ctmp = confra( ntrm, zinv, xx )
c
c *** downward recurrence for 'biga'
c *** ( ref. 1, eq. 22 )
if ( noabs ) then
c ** no-absorption case
rbiga( ntrm ) = dble( ctmp )
do 25 n = ntrm, 2, - 1
rbiga( n-1 ) = (n*rezinv)
$ - 1.0 / ( (n*rezinv) + rbiga( n ) )
25 continue
c
else
c ** absorptive case
cbiga( ntrm ) = ctmp
do 30 n = ntrm, 2, - 1
cbiga( n-1 ) = (n*zinv) - 1.0 / ( (n*zinv) + cbiga( n ) )
30 continue
c
end if
c
else
c *** upward recurrence for 'biga'
c *** ( ref. 1, eqs. 20-21 )
if ( noabs ) then
c ** no-absorption case
rtmp = dsin( mre*xx )
rbiga( 1 ) = - rezinv
$ + rtmp / ( rtmp*rezinv - dcos( mre*xx ) )
do 40 n = 2, ntrm
rbiga( n ) = - ( n*rezinv )
$ + 1.0 / ( ( n*rezinv ) - rbiga( n-1 ) )
40 continue
c
else
c ** absorptive case
ctmp = cdexp( - dcmplx(0.d0,2.d0) * cior * xx )
cbiga( 1 ) = - zinv + (1.-ctmp) / ( zinv * (1.-ctmp) -
$ dcmplx(0.d0,1.d0)*(1.+ctmp) )
do 50 n = 2, ntrm
cbiga( n ) = - (n*zinv) + 1.0 / ((n*zinv) - cbiga( n-1 ))
50 continue
end if
c
end if
c
return
end
complex*16 function confra( n, zinv, xx )
c
c compute bessel function ratio capital-a-sub-n from its
c continued fraction using lentz method ( ref. 1, pp. 17-20 )
c
c zinv = reciprocal of argument of capital-a
c
c i n t e r n a l v a r i a b l e s
c ------------------------------------
c
c cak term in continued fraction expansion of capital-a
c ( ref. 1, eq. 25 )
c capt factor used in lentz iteration for capital-a
c ( ref. 1, eq. 27 )
c cdenom denominator in -capt- ( ref. 1, eq. 28b )
c cnumer numerator in -capt- ( ref. 1, eq. 28a )
c cdtd product of two successive denominators of -capt-
c factors ( ref. 1, eq. 34c )
c cntn product of two successive numerators of -capt-
c factors ( ref. 1, eq. 34b )
c eps1 ill-conditioning criterion
c eps2 convergence criterion
c kk subscript k of -cak- ( ref. 1, eq. 25b )
c kount iteration counter ( used only to prevent runaway )
c maxit max. allowed no. of iterations
c mm + 1 and - 1, alternately
c
implicit none
integer n,maxit,mm,kk,kount
real*8 xx,eps1,eps2
complex*16 zinv
complex*16 cak, capt, cdenom, cdtd, cnumer, cntn
data eps1 / 1.e - 2 /, eps2 / 1.e - 8 /
data maxit / 10000 /
c
c *** ref. 1, eqs. 25a, 27
confra = ( n + 1 ) * zinv
mm = - 1
kk = 2 * n + 3
cak = ( mm * kk ) * zinv
cdenom = cak
cnumer = cdenom + 1.0 / confra
kount = 1
c
20 kount = kount + 1
if ( kount.gt.maxit )
$ call errmsg( 'confra--iteration failed to converge$', .true.)
c
c *** ref. 2, eq. 25b
mm = - mm
kk = kk + 2
cak = ( mm * kk ) * zinv
c *** ref. 2, eq. 32
if ( cdabs( cnumer/cak ).le.eps1
$ .or. cdabs( cdenom/cak ).le.eps1 ) then
c
c ** ill-conditioned case -- stride
c ** two terms instead of one
c
c *** ref. 2, eqs. 34
cntn = cak * cnumer + 1.0
cdtd = cak * cdenom + 1.0
confra = ( cntn / cdtd ) * confra
c *** ref. 2, eq. 25b
mm = - mm
kk = kk + 2
cak = ( mm * kk ) * zinv
c *** ref. 2, eqs. 35
cnumer = cak + cnumer / cntn
cdenom = cak + cdenom / cdtd
kount = kount + 1
go to 20
c
else
c ** well-conditioned case
c
c *** ref. 2, eqs. 26, 27
capt = cnumer / cdenom
confra = capt * confra
c ** check for convergence
c ** ( ref. 2, eq. 31 )
c
if ( dabs( dble(capt) - 1.0 ).ge.eps2
$ .or. dabs( dimag(capt) ) .ge.eps2 ) then
c
c *** ref. 2, eqs. 30a-b
cnumer = cak + 1.0 / cnumer
cdenom = cak + 1.0 / cdenom
go to 20
end if
end if
c
return
c
end
subroutine miprnt( prnt, xx, perfct, crefin, numang, xmu, 2
$ qext, qsca, gqsc, nmom, ipolzn, momdim,
$ calcmo, pmom, sforw, sback, tforw, tback,
$ s1, s2 )
c
c print scattering quantities of a single particle
c
implicit none
logical perfct, prnt(*), calcmo(*)
integer ipolzn, momdim, nmom, numang,i,m,j
real*8 gqsc, pmom( 0:momdim, * ), qext, qsca, xx, xmu(*)
real*8 fi1,fi2,fnorm
complex*16 crefin, sforw, sback, tforw(*), tback(*), s1(*), s2(*)
character*22 fmt
c
c
if ( perfct ) write ( 59, '(''1'',10x,a,1p,e11.4)' )
$ 'perfectly conducting case, size parameter =', xx
if ( .not.perfct ) write ( 59, '(''1'',10x,3(a,1p,e11.4))' )
$ 'refractive index: real ', dble(crefin),
$ ' imag ', dimag(crefin), ', size parameter =', xx
c
if ( prnt(1) .and. numang.gt.0 ) then
c
write ( 59, '(/,a)' )
$ ' cos(angle) ------- s1 --------- ------- s2 ---------'//
$ ' --- s1*conjg(s2) --- i1=s1**2 i2=s2**2 (i1+i2)/2'//
$ ' deg polzn'
do 10 i = 1, numang
fi1 = dble( s1(i) ) **2 + dimag( s1(i) )**2
fi2 = dble( s2(i) ) **2 + dimag( s2(i) )**2
write( 59, '( i4, f10.6, 1p,10e11.3 )' )
$ i, xmu(i), s1(i), s2(i), s1(i)*dconjg(s2(i)),
$ fi1, fi2, 0.5*(fi1+fi2), (fi2-fi1)/(fi2+fi1)
10 continue
c
end if
c
c
if ( prnt(2) ) then
c
write ( 59, '(/,a,9x,a,17x,a,17x,a,/,(0p,f7.2, 1p,6e12.3) )' )
$ ' angle', 's-sub-1', 't-sub-1', 't-sub-2',
$ 0.0, sforw, tforw(1), tforw(2),
$ 180., sback, tback(1), tback(2)
write ( 59, '(/,4(a,1p,e11.4))' )
$ ' efficiency factors, extinction:', qext,
$ ' scattering:', qsca,
$ ' absorption:', qext-qsca,
$ ' rad. pressure:', qext-gqsc
c
if ( nmom.gt.0 ) then
c
write( 59, '(/,a)' ) ' normalized moments of :'
if ( ipolzn.eq.0 ) write ( 59, '(''+'',27x,a)' ) 'phase fcn'
if ( ipolzn.gt.0 ) write ( 59, '(''+'',33x,a)' )
$ 'm1 m2 s21 d21'
if ( ipolzn.lt.0 ) write ( 59, '(''+'',33x,a)' )
$ 'r1 r2 r3 r4'
c
fnorm = 4. / ( xx**2 * qsca )
do 20 m = 0, nmom
write ( 59, '(a,i4)' ) ' moment no.', m
do 20 j = 1, 4
if( calcmo(j) ) then
write( fmt, 98 ) 24 + (j-1)*13
write ( 59,fmt ) fnorm * pmom(m,j)
end if
20 continue
end if
c
end if
c
return
c
98 format( '( ''+'', t', i2, ', 1p,e13.4 )' )
end
subroutine small1 ( xx, numang, xmu, qext, qsca, gqsc, sforw, 1
$ sback, s1, s2, tforw, tback, a, b )
c
c small-particle limit of mie quantities in totally reflecting
c limit ( mie series truncated after 2 terms )
c
c a,b first two mie coefficients, with numerator and
c denominator expanded in powers of -xx- ( a factor
c of xx**3 is missing but is restored before return
c to calling program ) ( ref. 2, p. 1508 )
c
implicit none
integer numang,j
real*8 gqsc, qext, qsca, xx, xmu(*)
real*8 twothr,fivthr,fivnin,sq,rtmp
complex*16 a( 2 ), b( 2 ), sforw, sback, s1(*), s2(*),
$ tforw(*), tback(*)
c
parameter ( twothr = 2./3., fivthr = 5./3., fivnin = 5./9. )
complex*16 ctmp
sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
c
c
a( 1 ) = dcmplx ( 0.d0, twothr * ( 1. - 0.2 * xx**2 ) )
$ / dcmplx ( 1.d0 - 0.5 * xx**2, twothr * xx**3 )
c
b( 1 ) = dcmplx ( 0.d0, - ( 1. - 0.1 * xx**2 ) / 3. )
$ / dcmplx ( 1.d0 + 0.5 * xx**2, - xx**3 / 3. )
c
a( 2 ) = dcmplx ( 0.d0, xx**2 / 30. )
b( 2 ) = dcmplx ( 0.d0, - xx**2 / 45. )
c
qsca = 6. * xx**4 * ( sq( a(1) ) + sq( b(1) )
$ + fivthr * ( sq( a(2) ) + sq( b(2) ) ) )
qext = qsca
gqsc = 6. * xx**4 * dble( a(1) * dconjg( a(2) + b(1) )
$ + ( b(1) + fivnin * a(2) ) * dconjg( b(2) ) )
c
rtmp = 1.5 * xx**3
sforw = rtmp * ( a(1) + b(1) + fivthr * ( a(2) + b(2) ) )
sback = rtmp * ( a(1) - b(1) - fivthr * ( a(2) - b(2) ) )
tforw( 1 ) = rtmp * ( b(1) + fivthr * ( 2.*b(2) - a(2) ) )
tforw( 2 ) = rtmp * ( a(1) + fivthr * ( 2.*a(2) - b(2) ) )
tback( 1 ) = rtmp * ( b(1) - fivthr * ( 2.*b(2) + a(2) ) )
tback( 2 ) = rtmp * ( a(1) - fivthr * ( 2.*a(2) + b(2) ) )
c
do 10 j = 1, numang
s1( j ) = rtmp * ( a(1) + b(1) * xmu(j) + fivthr *
$ ( a(2) * xmu(j) + b(2) * ( 2.*xmu(j)**2 - 1. )) )
s2( j ) = rtmp * ( b(1) + a(1) * xmu(j) + fivthr *
$ ( b(2) * xmu(j) + a(2) * ( 2.*xmu(j)**2 - 1. )) )
10 continue
c ** recover actual mie coefficients
a( 1 ) = xx**3 * a( 1 )
a( 2 ) = xx**3 * a( 2 )
b( 1 ) = xx**3 * b( 1 )
b( 2 ) = xx**3 * b( 2 )
c
return
end
subroutine small2 ( xx, cior, calcqe, numang, xmu, qext, qsca, 1
$ gqsc, sforw, sback, s1, s2, tforw, tback,
$ a, b )
c
c small-particle limit of mie quantities for general refractive
c index ( mie series truncated after 2 terms )
c
c a,b first two mie coefficients, with numerator and
c denominator expanded in powers of -xx- ( a factor
c of xx**3 is missing but is restored before return
c to calling program ) ( ref. 2, p. 1508 )
c
c ciorsq square of refractive index
c
implicit none
logical calcqe
integer numang,j
real*8 gqsc, qext, qsca, xx, xmu(*)
real*8 twothr,fivthr,sq,rtmp
complex*16 a( 2 ), b( 2 ), cior, sforw, sback, s1(*), s2(*),
$ tforw(*), tback(*)
c
parameter ( twothr = 2./3., fivthr = 5./3. )
complex*16 ctmp, ciorsq
sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
c
c
ciorsq = cior**2
ctmp = dcmplx( 0.d0, twothr ) * ( ciorsq - 1.0 )
a(1) = ctmp * ( 1.0 - 0.1 * xx**2 + (ciorsq/350. + 1./280.)*xx**4)
$ / ( ciorsq + 2.0 + ( 1.0 - 0.7 * ciorsq ) * xx**2
$ - ( ciorsq**2/175. - 0.275 * ciorsq + 0.25 ) * xx**4
$ + xx**3 * ctmp * ( 1.0 - 0.1 * xx**2 ) )
c
b(1) = (xx**2/30.) * ctmp * ( 1.0 + (ciorsq/35. - 1./14.) *xx**2 )
$ / ( 1.0 - ( ciorsq/15. - 1./6. ) * xx**2 )
c
a(2) = ( 0.1 * xx**2 ) * ctmp * ( 1.0 - xx**2 / 14. )
$ / ( 2. * ciorsq + 3. - ( ciorsq/7. - 0.5 ) * xx**2 )
c
qsca = 6. * xx**4 * ( sq(a(1)) + sq(b(1)) + fivthr * sq(a(2)) )
gqsc = 6. * xx**4 * dble( a(1) * dconjg( a(2) + b(1) ) )
qext = qsca
if ( calcqe ) qext = 6. * xx * dble( a(1) + b(1) + fivthr * a(2) )
c
rtmp = 1.5 * xx**3
sforw = rtmp * ( a(1) + b(1) + fivthr * a(2) )
sback = rtmp * ( a(1) - b(1) - fivthr * a(2) )
tforw( 1 ) = rtmp * ( b(1) - fivthr * a(2) )
tforw( 2 ) = rtmp * ( a(1) + 2. * fivthr * a(2) )
tback( 1 ) = tforw( 1 )
tback( 2 ) = rtmp * ( a(1) - 2. * fivthr * a(2) )
c
do 10 j = 1, numang
s1( j ) = rtmp * ( a(1) + ( b(1) + fivthr * a(2) ) * xmu(j) )
s2( j ) = rtmp * ( b(1) + a(1) * xmu(j) + fivthr * a(2)
$ * ( 2. * xmu(j)**2 - 1. ) )
10 continue
c ** recover actual mie coefficients
a( 1 ) = xx**3 * a( 1 )
a( 2 ) = xx**3 * a( 2 )
b( 1 ) = xx**3 * b( 1 )
b( 2 ) = ( 0., 0. )
c
return
end
subroutine testmi ( qext, qsca, gqsc, sforw, sback, s1, s2, 1
$ tforw, tback, pmom, momdim, ok )
c
c compare mie code test case results with correct answers
c and return ok=false if even one result is inaccurate.
c
c the test case is : mie size parameter = 10
c refractive index = 1.5 - 0.1 i
c scattering angle = 140 degrees
c 1 sekera moment
c
c results for this case may be found among the test cases
c at the end of reference (1).
c
c *** note *** when running on some computers, esp. in single
c precision, the 'accur' criterion below may have to be relaxed.
c however, if 'accur' must be set larger than 10**-3 for some
c size parameters, your computer is probably not accurate
c enough to do mie computations for those size parameters.
c
implicit none
integer momdim,m,n
real*8 qext, qsca, gqsc, pmom( 0:momdim, * )
complex*16 sforw, sback, s1(*), s2(*), tforw(*), tback(*)
logical ok, wrong
c
real*8 accur, testqe, testqs, testgq, testpm( 0:1 )
complex*16 testsf, testsb,tests1,tests2,testtf(2), testtb(2)
data testqe / 2.459791 /, testqs / 1.235144 /,
$ testgq / 1.139235 /, testsf / ( 61.49476, -3.177994 ) /,
$ testsb / ( 1.493434, 0.2963657 ) /,
$ tests1 / ( -0.1548380, -1.128972) /,
$ tests2 / ( 0.05669755, 0.5425681) /,
$ testtf / ( 12.95238, -136.6436 ), ( 48.54238, 133.4656 ) /,
$ testtb / ( 41.88414, -15.57833 ), ( 43.37758, -15.28196 )/,
$ testpm / 227.1975, 183.6898 /
real*8 calc,exact
c data accur / 1.e-5 /
data accur / 1.e-4 /
wrong( calc, exact ) = dabs( (calc - exact) / exact ) .gt. accur
c
c
ok = .true.
if ( wrong( qext,testqe ) )
$ call tstbad( 'qext', abs((qext - testqe) / testqe), ok )
if ( wrong( qsca,testqs ) )
$ call tstbad( 'qsca', abs((qsca - testqs) / testqs), ok )
if ( wrong( gqsc,testgq ) )
$ call tstbad( 'gqsc', abs((gqsc - testgq) / testgq), ok )
c
if ( wrong( dble(sforw), dble(testsf) ) .or.
$ wrong( dimag(sforw), dimag(testsf) ) )
$ call tstbad( 'sforw', cdabs((sforw - testsf) / testsf), ok )
c
if ( wrong( dble(sback), dble(testsb) ) .or.
$ wrong( dimag(sback), dimag(testsb) ) )
$ call tstbad( 'sback', cdabs((sback - testsb) / testsb), ok )
c
if ( wrong( dble(s1(1)), dble(tests1) ) .or.
$ wrong( dimag(s1(1)), dimag(tests1) ) )
$ call tstbad( 's1', cdabs((s1(1) - tests1) / tests1), ok )
c
if ( wrong( dble(s2(1)), dble(tests2) ) .or.
$ wrong( dimag(s2(1)), dimag(tests2) ) )
$ call tstbad( 's2', cdabs((s2(1) - tests2) / tests2), ok )
c
do 20 n = 1, 2
if ( wrong( dble(tforw(n)), dble(testtf(n)) ) .or.
$ wrong( dimag(tforw(n)), dimag(testtf(n)) ) )
$ call tstbad( 'tforw', cdabs( (tforw(n) - testtf(n)) /
$ testtf(n) ), ok )
if ( wrong( dble(tback(n)), dble(testtb(n)) ) .or.
$ wrong( dimag(tback(n)), dimag(testtb(n)) ) )
$ call tstbad( 'tback', cdabs( (tback(n) - testtb(n)) /
$ testtb(n) ), ok )
20 continue
c
do 30 m = 0, 1
if ( wrong( pmom(m,1), testpm(m) ) )
$ call tstbad( 'pmom', dabs( (pmom(m,1)-testpm(m)) /
$ testpm(m) ), ok )
30 continue
c
return
c
end
subroutine errmsg( messag, fatal ) 9
c
c print out a warning or error message; abort if error
c
implicit none
logical fatal, once
character*(*) messag
integer maxmsg, nummsg
save maxmsg, nummsg, once
data nummsg / 0 /, maxmsg / 100 /, once / .false. /
c
c
if ( fatal ) then
write ( *, '(2a)' ) ' ******* error >>>>>> ', messag
stop
end if
c
nummsg = nummsg + 1
if ( nummsg.gt.maxmsg ) then
if ( .not.once ) write ( *,99 )
once = .true.
else
write ( *, '(2a)' ) ' ******* warning >>>>>> ', messag
endif
c
return
c
99 format( ///,' >>>>>> too many warning messages -- ',
$ 'they will no longer be printed <<<<<<<', /// )
end
subroutine wrtbad ( varnam, erflag ) 9
c
c write names of erroneous variables
c
c input : varnam = name of erroneous variable to be written
c ( character, any length )
c
c output : erflag = logical flag, set true by this routine
c ----------------------------------------------------------------------
implicit none
character*(*) varnam
logical erflag
integer maxmsg, nummsg
save nummsg, maxmsg
data nummsg / 0 /, maxmsg / 50 /
c
c
nummsg = nummsg + 1
write ( *, '(3a)' ) ' **** input variable ', varnam,
$ ' in error ****'
erflag = .true.
if ( nummsg.eq.maxmsg )
$ call errmsg ( 'too many input errors. aborting...$', .true. )
return
c
end
subroutine tstbad( varnam, relerr, ok ) 10
c
c write name (-varnam-) of variable failing self-test and its
c percent error from the correct value. return ok = false.
c
implicit none
character*(*) varnam
logical ok
real*8 relerr
c
c
ok = .false.
write( *, '(/,3a,1p,e11.2,a)' )
$ ' output variable ', varnam,' differed by', 100.*relerr,
$ ' per cent from correct value. self-test failed.'
return
c
end