c
c THE MILNE PROBLEM
c =================================
c
c The following codes constitute a package of subroutines in FORTRAN-77
c for finding the radiation field in a homogeneous optically semi-infinite
c plane-parallel atmosphere where the scattering obeys the
c Rayleigh-Cabannes law and the sources of radiation are infinitely
c deep in the atmosphere.
c
c The theoretical background is given in papers:
c
c 1. Viik T., RAYLEIGH-CABANNES SCATTERING IN PLANETARY ATMOSPHERES
c III. THE MILNE PROBLEM IN CONSERVATIVE ATMOSPHERES,
c Earth, Moon, and Planets, vol. 49, 163-175, 1990.
c
c 2. Viik T., RADIATION FIELD IN RAYLEIGH-CABANNES SCATTERING
c ATMOSPHERE: THE NONCONSERVATIVE MILNE PROBLEM,
c J. Quant. Spectr. Radiative Transfer, vol. 66,
c 581-590, 2000.
c
c
c The main input parameters are:
c alamb - albedo of single scattering
c qs1 - parameter of depolarization
c (qs1=2*(1-\rho_n)/(2+\rho_n), where
c \rho_n is the depolarization factor)
c nq - the order of the Gauss quadrature
c
c For checking purposes I have included the code
c CHECK_MILNE.FOR
c together with a file of results
c RESULT_MILNE.DAT.
c
c
c
c GENERAL CODES
c
c
c
c GAULEG - finds the array of points x(1...n) and the array of weights
c w(1...n) for numerical integration in interval (x1,x2).
c Source - Numerical Recipes.
c
c LINSYS - solves a set of linear algebraic equations while the array
c of coefficients is a(n*n) and the right-hand vector is b(n).
c The solution is stored in the array b.
c Source - IBM software.
c
c APPRCHAF - evaluates the approximate characteristic function at x
c for Rayleigh-Cabannes phase matrix (2 by 2).
c
c
c ZEROPOL - evaluates the zeros of the approximate
c characteristic equation for Rayleigh-Cabannes
c phase matrix (2 by 2) and stores them in arrays
c sbr(nq) and sbl(nq).
c
c ZERO1 - finds zero of function f if found in the
c interval ax,bx.
c Source - G.E.Forsythe, M.A.Malcolm, C.B.Moler
c
c ZERO2 - finds zero of function f if found in the
c interval ax,bx.
c Source - G.E.Forsythe, M.A.Malcolm, C.B.Moler
c
c CHARFUNCT - finds the characteristic function
c
c CHAREQ - evaluates the approximate characteristic equation
c as a function of x
c
c ZERONAT - evaluates the zeros of the approximate characteristic
c equation and stores them in array s
c
c
c DK024 - finds for Rayleigh-Cabannes scattering the values of
c DA_0, DA_2, DA_4, DB_0, DB_2 and DB_4
c
c
c CODES FOR CONSERVATIVE SCATTERING
c
c
c ALBEGAK - finds for Rayleigh-Cabannes scattering the values of
c functions z(k) and v(k)
c
c MILNCO - finds the coefficients a(i) and c(i) for 2*2
c Rayleigh-Cabannes scattering
c
c MINTG - finds intensities at any level for Rayleigh-Cabannes
c scattering in homogeneous plane-parallel conservative
c atmosphere at the Gaussian points uu(i), i=1,...,nq
c
c MINT - finds the radiation field for a 2*2 Rayleigh-Cabannes
c scattering atmosphere at the angle points
c us(i) (i=1,...,ni)
c
c MSSQQ - finds the source functions s_i and s_q together with
c q_i and q_q functions for the 2*2 Rayleigh-Cabannes
c scattering atmosphere
c
c
c CODES FOR NONCONSERVATIVE SCATTERING
c
c
c
c WURZEL - finds the solution of the characteristic equation
c for vector transfer
c
c RFX - finds the left hand side of the characteristic
c equation for the Rayleigh-Cabannes scattering
c
c RTBISX - using bisection, finds the root of a function func known
c to lie between x1 and x2.
c Source - Numerical Recipes
c
c ZAVA - finds for Rayleigh scattering the values of
c functions z(k) and v(k)
c
c COEFMI - finds the coefficients a(i) and c(i)
c for the Milne problem
c
c MINTGN - for the non-conservative Milne problem it finds
c the intensities at any level for 2*2
c Rayleigh-Cabannes scattering atmosphere at
c Gaussian points uu(i), i=1,...,nq.
c
c MNCINT - for the non-conservative Milne problem it finds
c the intensities at any level for 2*2
c Rayleigh-Cabannes scattering atmosphere at
c arbitrary point us(i), i=1,...,ni.
c
c MNCQQ - for the non-conservative Milne problem it finds
c the q_i and q_q functions for the 2*2
c Rayleigh-Cabannes scattering atmosphere
c
c MNCSOU - for the non-conservative Milne problem it finds
c the s_l and s_r functions together with s_i and s_q
c functions for the 2*2 Rayleigh-Cabannes scattering
c atmosphere
c
c
c
c GENERAL CODES
c
c
subroutine gauleg(x1,x2,x,w,n)
implicit real*8(a-h,o-z)
c given the lower and upper limits of integration x1 and x2
c and given n, this code returns arrays x and w of length n
c containing the abscissas and weigths of the Gauss-
c Legendre n-point quadrature formula.
parameter (eps=.11d-15,pi=3.141592653589793238d0)
dimension x(n),w(n)
m=(n+1)/2
itmax=40
xm=.5d0*(x2+x1)
xl=.5d0*(x2-x1)
do 12 i=1,m
z =dcos(pi*(dfloat(i)-.25d0)/(dfloat(n)+.5d0))
it=0
1 p1=1.d0
p2=0.d0
do 11 j=1,n
p3=p2
p2=p1
p1=((2.d0*dfloat(j)-1.d0)*z*p2-(dfloat(j)-1.d0)*p3)/
& dfloat(j)
11 continue
pp=dfloat(n)*(z*p1-p2)/(z*z-1.d0)
z1=z
z=z1-p1/pp
it=it+1
if (dabs(z-z1).gt.eps.and.it.lt.itmax) go to 1
x(i)=xm-xl*z
x(n+1-i)=xm+xl*z
w(i)=2.d0*xl/((1.d0-z*z)*pp*pp)
w(n+1-i)=w(i)
12 continue
return
end
subroutine linsys(a,b,n,ks)
implicit real*8(a-h,o-z)
c solution of linear system ax=b
c a is stored as follows: 1.column,2.column etc.
c b is rhs vector of order n
c n=number of unknowns (and equations)
c solution is stored in b
c ks=output index
c =0_ normal solution
c =1_ singular matrix
dimension a(1),b(1)
c 1
tol=0d0
ks=0
jj=-n
do 65 j=1,n
jy=j+1
jj=jj+n+1
biga=0d0
it=jj-j
do 30 i=j,n
c 2
ij=it+i
if (dabs(biga)-dabs(a(ij))) 20,30,30
20 biga=a(ij)
imax=i
30 continue
c 3
if (dabs(biga)-tol) 35,35,40
35 ks=1
return
c 4
40 ii=j+n*(j-2)
it=imax-j
do 50 k=j,n
ii=ii+n
i2=ii+it
save=a(ii)
a(ii)=a(i2)
a(i2)=save
c 5
50 a(ii)=a(ii)/biga
save=b(imax)
b(imax)=b(j)
b(j)=save/biga
c 6
if(j-n) 55,70,55
55 iqs=n*(j-1)
do 65 ix=jy,n
ixj=iqs+ix
it=j-ix
do 60 jx=jy,n
ixjx=n*(jx-1)+ix
jjx=ixjx+it
60 a(ixjx)=a(ixjx)-(a(ixj)*a(jjx))
65 b(ix)=b(ix)-(b(j)*a(ixj))
c 7
70 ny=n-1
it=n*n
do 80 j=1,ny
ia=it-j
ib=n-j
ic=n
do 80 k=1,j
b(ib)=b(ib)-a(ia)*b(ic)
ia=ia-n
80 ic=ic-1
return
end
subroutine zeropol(sbr,sbl,alamb,qs1,uu,cc,nq)
implicit real*8(a-h,o-z)
parameter (nmax=400,sioux=1d-10,nqmax=200)
c this code evaluates the zeros of the approximate
c characteristic equation for Rayleigh-Cabannes phase matrix (2 by 2)
c and stores them in arrays sbr(nq) and sbl(nq).
c Viik 87 12 06
logical lg1,lg2,lg3,lg4,m1,m2
dimension uu(nq),cc(nq),sbr(nq),sbl(nq),s(nqmax)
external apprchaf
data tol/1.d-14/
tolex=0.23d-15
itmax=50
stol=1.d0-tol
m1=alamb.ge.stol
m2=qs1.ge.stol
if (m1.and.m2) then
do 20 liw=1,2
miw=liw+8
call zeronat(s,alamb,uu,cc,nq,miw)
do 30 i=1,nq
if (liw.eq.1) sbr(i)=s(i)
if (liw.eq.2) sbl(i)=s(i)
30 continue
20 continue
return
endif
do 100 ij=1,nq
ji=nq+1-ij
u=1.d0/uu(ji)
it=0
ax=sioux
bx=u-sioux
if (ij.ne.1) ax=1.d0/uu(ji+1)+sioux
if (ij.eq.1.and.m1) go to 60
40 if (it.gt.itmax) go to 100
cx=(ax+bx)/2.d0
fax=apprchaf(ax,alamb,qs1,uu,cc,nq)
fcx=apprchaf(cx,alamb,qs1,uu,cc,nq)
lg1=fax.gt.0.d0
lg2=fcx.gt.0.d0
if (lg1.and.lg2.or..not.lg1.and..not.lg2) go to 50
sbl(ij)=zero1(ax,cx,apprchaf,tolex,alamb,qs1,uu,cc,nq)
sbr(ij)=zero1(cx,bx,apprchaf,tolex,alamb,qs1,uu,cc,nq)
go to 100
50 dx=(cx-ax)/dfloat(nmax)
dax=ax+dx
fdax=apprchaf(dax,alamb,qs1,uu,cc,nq)
dcx=cx-dx
fdcx=apprchaf(dcx,alamb,qs1,uu,cc,nq)
da=(fdax-fax)/dx
dc=(fcx-fdcx)/dx
lg3=da.gt.0.d0
lg4=dc.gt.0.d0
if (lg3.and.lg4.or..not.lg3.and..not.lg4) go to 70
bx=cx
it=it+1
go to 40
70 ax=cx
it=it+1
go to 40
60 ax=.5d0*(ax+bx)
sbl(1)=0.d0
sbr(1)=zero1(ax,bx,apprchaf,tolex,alamb,qs1,uu,cc,nq)
100 continue
return
end
function apprchaf(x,alamb,qs1,uu,cc,nq)
implicit real*8(a-h,o-z)
c apprchaf evaluates the approximate
c characteristic function at x for Rayleigh phase matrix (2 by 2).
c Viik 87 12 06
dimension uu(nq),cc(nq)
data f1/1.333333333333333d0/,f2/2.666666666666667d0/,
& f3/.888888888888889d0/
d0=0.d0
d2=0.d0
d4=0.d0
x2=x*x
do 10 i=1,nq
c=cc(i)
u=uu(i)
u2=u*u
u4=u2*u2
h=1.d0-x2*u2
20 ch=c/h
d0=d0+ch
d2=d2+ch*u2
d4=d4+ch*u4
10 continue
d0=2.d0*d0
d2=2.d0*d2
d4=2.d0*d4
z=alamb*(d0-d2)
qs2=1.d0-qs1
v=(z-f1)*(z*qs1-f2)
v=v-4.d0*alamb*qs1*(1.d0-alamb)*d4
v=v+f1*alamb*d2*(qs1*(3.d0-alamb)-2.d0)
apprchaf=v+f3*alamb*qs2*d0
return
end
function zero1(ax,bx,f,tol,alamb,qs1,uu,cc,nq)
implicit real*8(a-h,o-z)
c zero1 determines zero of function f if found in the
c interval ax,bx
c nq is the order of the Gaussian quadrature
c uu(nq) is the array of Gaussian points
c cc(nq) is the array of Gaussian weights
c tol - desired length of the interval
c of uncertainty for the zero
c G.E.Forsythe, M.A.Malcolm, C.B.Moler
c Viik sept 1980
dimension uu(nq),cc(nq)
data eps/.222d-15/
a=ax
b=bx
fa=f(a,alamb,qs1,uu,cc,nq)
fb=f(b,alamb,qs1,uu,cc,nq)
20 c=a
fc=fa
d=b-a
e=d
30 if (dabs(fc).ge.dabs(fb)) go to 40
a=b
b=c
c=a
fa=fb
fb=fc
fc=fa
40 tol1=2.d0*eps*dabs(b)+0.5d0*tol
xm=0.5d0*(c-b)
if (dabs(xm).le.tol1) go to 90
if (fb.eq.0.d0) go to 90
if (dabs(e).lt.tol1) go to 70
if (dabs(fa).le.dabs(fb)) go to 70
if (a.ne.c) go to 50
s=fb/fa
p=2.d0*xm*s
q=1.d0-s
go to 60
50 q=fa/fc
r=fb/fc
s=fb/fa
p=s*(2.d0*xm*q*(q-r)-(b-a)*(r-1.d0))
q=(q-1.d0)*(r-1.d0)*(s-1.d0)
60 if (p.gt.0.d0) q=-q
p=dabs(p)
if ((2.d0*p).ge.(3.d0*xm*q-dabs(tol1*q))) go to 70
if (p.ge.dabs(0.5d0*e*q)) go to 70
e=d
d=p/q
go to 80
70 d=xm
e=d
80 a=b
fa=fb
if (dabs(d).gt.tol1) b=b+d
if (dabs(d).le.tol1) b=b+dsign(tol1,xm)
fb=f(b,alamb,qs1,uu,cc,nq)
if ((fb*(fc/dabs(fc))).gt.0.d0) go to 20
go to 30
90 zero1=b
return
end
function zero2(ax,bx,f,tol,alamb,uu,cc,nq,miw)
implicit real*8(a-h,o-z)
c zero2 determines zero of function f if found in the
c interval ax,bx
c nq is the order of the Gaussian quadrature
c uu(nq) is the array of Gaussian points
c cc(nq) is the array of Gaussian weights
c tol - desired length of the interval
c of uncertainty for the zero
c miw chooses the proper characteristic function
c G.E.Forsythe, M.A.Malcolm, C.B.Moler
c Viik sept 1980
dimension uu(nq),cc(nq)
data eps/.222d-15/
a=ax
b=bx
fa=f(a,alamb,uu,cc,nq,miw)
fb=f(b,alamb,uu,cc,nq,miw)
20 c=a
fc=fa
d=b-a
e=d
30 if (dabs(fc).ge.dabs(fb)) go to 40
a=b
b=c
c=a
fa=fb
fb=fc
fc=fa
40 tol1=2.d0*eps*dabs(b)+0.5d0*tol
xm=0.5d0*(c-b)
if (dabs(xm).le.tol1) go to 90
if (fb.eq.0.d0) go to 90
if (dabs(e).lt.tol1) go to 70
if (dabs(fa).le.dabs(fb)) go to 70
if (a.ne.c) go to 50
s=fb/fa
p=2.d0*xm*s
q=1.d0-s
go to 60
50 q=fa/fc
r=fb/fc
s=fb/fa
p=s*(2.d0*xm*q*(q-r)-(b-a)*(r-1.d0))
q=(q-1.d0)*(r-1.d0)*(s-1.d0)
60 if (p.gt.0.d0) q=-q
p=dabs(p)
if ((2.d0*p).ge.(3.d0*xm*q-dabs(tol1*q))) go to 70
if (p.ge.dabs(0.5d0*e*q)) go to 70
e=d
d=p/q
go to 80
70 d=xm
e=d
80 a=b
fa=fb
if (dabs(d).gt.tol1) b=b+d
if (dabs(d).le.tol1) b=b+dsign(tol1,xm)
fb=f(b,alamb,uu,cc,nq,miw)
if ((fb*(fc/dabs(fc))).gt.0.d0) go to 20
go to 30
90 zero2=b
return
end
function charfunct(iw,u,alamb)
implicit real*8(a-h,o-z)
c charfunct is the characteristic function
c
c iw=0 - isotr. scat.
c
c Rayleigh scattering (with Mullikin indeces)
c
c iw=6 - Mullikin index 1
c iw=7 - Mullikin index 2
c iw=8 - Mullikin index 3
c iw=9 - Mullikin index 4
c iw=10 - Mullikin index 5
c u - angle variable,
c alamb - albedo of single scattering
c
c This code is used by CHAREQ
c
c Viik 84 02 20
u2=u*u
uv=1.d0-u2
al=1.d0-alamb
if (iw.eq.0) charfunct=.5d0*alamb
if (iw.eq.6) charfunct=0.375d0*alamb*uv*(1.d0+2.d0*u2)
if (iw.eq.7) charfunct=0.1875d0*alamb*(1.d0+u2)*(1.d0+u2)
if (iw.eq.8) charfunct=0.75d0*alamb*u2
if (iw.eq.9) charfunct=0.375d0*alamb*uv
if (iw.eq.10) charfunct=0.75d0*alamb*uv
return
end
function chareq(x,alamb,uu,cc,nq,miw)
implicit real*8(a-h,o-z)
c chareq evaluates the approximate characteristic equation
c as a function of x
c x - argument on the beta-axis
c This code is used by ZERONAT
dimension uu(nq),cc(nq)
external charfunct
iw=miw
sum=0.d0
do 100 i=1,nq
u=uu(i)
u2=u*u
w=charfunct(iw,u,alamb)
sx=1.d0/(1.d0+x*u)+1.d0/(1.d0-x*u)
sum=sum+cc(i)*w*sx
100 continue
chareq=1.d0-sum
return
end
subroutine zeronat(s,alamb,uu,cc,nq,miw)
implicit real*8(a-h,o-z)
c zeronat evaluates the zeros of the
c approximate characteristic equation
c and stores them in array s
logical lg1,lg2
dimension uu(nq),cc(nq),s(nq)
external chareq
data siu/1.d-8/,tol/.222d-15/
c siu - value by which the starting value
c of s is shifted from the
c asymptotic value in order to avoid
c overflow
c tol = computer-dependent constant,
c it is the smallest positive number such
c that (1.0+tol.gt.1.0)
tolex=0.d0
itmax=50
stol=1.d0-tol
niw=0
lg1=miw.eq.0.or.miw.eq.1
lg2=miw.eq.3.or.miw.eq.10
if (lg1.or.lg2) niw=1
do 100 ij=1,nq
if (niw.eq.1.and.alamb.gt.stol.and.ij.eq.1) go to 100
ji=nq+1-ij
u=1.d0/uu(ji)
ax=siu
bx=u-siu
if (ij.ne.1) ax=1.d0/uu(ji+1)+siu
s(ij)=zero2(ax,bx,chareq,tolex,alamb,uu,cc,nq,miw)
c zero2 = function subroutine for determining the zeros,
C e.g. - G.E.Forsythe, M.A.Malcolm, C.B.Moler
C COMPUTER METHODS FOR MATHEMATICAL CALCULATIONS.
C ENGLEWOOD CLIFFS, N.J.07632, PRENTICE-HALL, INC., 1977
100 continue
if (niw.eq.1.and.alamb.gt.stol) s(1)=0.d0
return
end
subroutine dk024(dk0,dk2,dk4,da0,da2,da4,db0,db2,db4,x,
& sbr,sbl,uu,cc,nq)
implicit real*8(a-h,o-z)
c this code finds for Rayleigh scattering
c the values of da0,da2,da4,db0,db2 and db4
c
c The Milne problem
c
c Viik 89 12 09
dimension da0(nq),da2(nq),da4(nq),db0(nq),db2(nq),db4(nq),
& cc(nq),uu(nq),sbr(nq),sbl(nq)
do 10 ix=1,2
do 30 i=1,nq
be=sbr(i)
if (ix.eq.2) be=sbl(i)
d0=0.d0
d2=0.d0
d4=0.d0
do 40 j=1,nq
u2=uu(j)*uu(j)
u4=u2*u2
be2=u2*be*be
cj=cc(j)/(1.d0-be2)
d0=d0+cj
d2=d2+cj*u2
40 d4=d4+cj*u4
d0=2.d0*d0
d2=2.d0*d2
d4=2.d0*d4
if (ix.ne.1) go to 50
da0(i)=d0
da2(i)=d2
da4(i)=d4
go to 30
50 db0(i)=d0
db2(i)=d2
db4(i)=d4
30 continue
10 continue
dk0=0.d0
dk2=0.d0
dk4=0.d0
do 20 j=1,nq
u2=uu(j)*uu(j)
u4=u2*u2
be2=u2*x*x
cj=cc(j)/(1.d0-be2)
dk0=dk0+cj
dk2=dk2+cj*u2
20 dk4=dk4+cj*u4
dk0=2.d0*dk0
dk2=2.d0*dk2
dk4=2.d0*dk4
return
end
c
c
c CONSERVATIVE CASE
c
c
subroutine albegak(zal,val,zbe,vbe,da0,da2,da4,
& db0,db2,db4,qs1,nq)
implicit real*8(a-h,o-z)
c this code finds for Rayleigh scattering
c the values of functions z(k) and v(k)
c
c The Milne problem
c
c Viik 89 12 09
dimension zal(nq),val(nq),zbe(nq),vbe(nq),
& da0(nq),da2(nq),da4(nq),db0(nq),db2(nq),db4(nq)
if (qs1.gt.(1.d0-1.d-12)) then
do 20 i=1,nq
zal(i)=0.d0
val(i)=0.d0
zbe(i)=0.d0
vbe(i)=0.d0
20 continue
return
endif
qs2=1.d0-qs1
aks=.375d0*qs1
dze=.25d0*qs2
do 10 ix=1,2
do 30 i=1,nq
if (ix.eq.1) then
d0=da0(i)
d2=da2(i)
d4=da4(i)
else
d0=db0(i)
d2=db2(i)
d4=db4(i)
endif
a11=1.d0-2.d0*aks*(d0-d2)-dze*d0
a12=-dze*d0
a21=-aks*d2-dze*d0
a22=1.d0-d0*(aks+dze)
b1=2.d0*aks*(d2-d4)+dze*d2
b2=aks*d4+dze*d2
del=a11*a22-a12*a21
z=(b1*a22-b2*a12)/del
v=(b2*a11-b1*a21)/del
if (ix.ne.1) go to 50
zal(i)=z
val(i)=v
go to 30
50 zbe(i)=z
vbe(i)=v
30 continue
10 continue
return
end
subroutine milnco(bm,zal,val,zbe,vbe,am,qs1,sbr,sbl,uu,nq)
implicit real*8(a-h,o-z)
c this code finds the coefficients a(i) and c(i)
c for 2*2 Rayleigh-Cabannes' scattering
c
c The Milne problem
c
c The scattering is conservative!!!!!!!!!
c
c Viik 89 12 09
c
c
logical ll
dimension jmin(4),jmax(4),imin(4),imax(4),zal(1),val(1),
& zbe(1),vbe(1),am(1),bm(1),sbr(nq),sbl(nq),uu(nq)
ll=qs1.le.(1.d0-1.d-12)
n2=2*nq
nl=nq
ij=0
do 131 i=1,2
do 131 j=1,2
ij=ij+1
imin(ij)=(i-1)*nq+1
imax(ij)=i*nq
jmin(ij)=(j-1)*nq+1
131 jmax(ij)=j*nq
do 132 ki=1,4
j1=jmin(ki)
j2=jmax(ki)
i1=imin(ki)
i2=imax(ki)
do 133 i=i1,i2
ie=i-i1+1
u=uu(ie)
u2=u*u
uv=1.d0-u2
do 133 j=j1,j2
je=j-j1+1
ij=i+n2*(j-1)
u5=1.d0-u*sbl(je)
br=sbr(je)
u4=1.d0-u*br
if (ll) then
if (ki.eq.1) am(ij)=(u2+zbe(je))/u5
if (ki.eq.2) am(ij)=(u2+zal(je))/u4
if (ki.eq.3) am(ij)=vbe(je)/u5
if (ki.eq.4) am(ij)=val(je)/u4
else if (.not.ll) then
if (ki.eq.1) am(ij)=uv/u5
if (ki.eq.2) am(ij)=1.d0+br*u
if (ki.eq.3) am(ij)=0.d0
if (ki.eq.4) am(ij)=(1.d0-br*br)/u4
endif
if (ki.eq.1.and.je.eq.1) am(ij)=1.d0
if (ki.eq.3.and.je.eq.1) am(ij)=1.d0
133 continue
132 continue
do 134 i=1,n2
if (i.le.nq) bm(i)=uu(i)
if (i.gt.nq) bm(i)=uu(i-nq)
134 continue
call linsys(am,bm,n2,kes)
return
end
subroutine mintg(aiul,aidl,aiur,aidr,tau,qs1,bm,zal,val,
& zbe,vbe,sbr,sbl,uu,nq)
implicit real*8(a-h,o-z)
parameter (nqmax=200)
c this code finds intensities at any level for
c 2*2 Rayleigh-Cabannes scattering in homogeneous
c plane-parallel conservative atmosphere
c at the Gaussian points
c
c The Milne problem
c
c Viik 89 12 09
logical ll
dimension bm(1),aiul(1),aidl(1),sbr(nq),sbl(nq),uu(nq),
& aiur(1),aidr(1),zal(1),val(1),zbe(1),vbe(1),
& e51(nqmax),e41(nqmax)
ll=qs1.le.(1.d0-1.d-12)
do 130 k=1,nq
bl=sbl(k)
br=sbr(k)
e51(k)=dexp(-bl*tau)
e41(k)=dexp(-br*tau)
130 continue
do 100 i=1,nq
v=uu(i)
sp1=0.d0
sm1=0.d0
sp2=0.d0
sm2=0.d0
sp3=0.d0
sm3=0.d0
sp4=0.d0
sm4=0.d0
do 110 j=2,nq
bl=sbl(j)
e1=e51(j)
vp=1.d0+v*bl
vm=1.d0-v*bl
vzb=v*v+zbe(j)
vv1=1.d0-v*v
if (ll) then
cm1=vzb/vm
cp1=vzb/vp
cm3=vbe(j)/vm
cp3=vbe(j)/vp
else
cm1=vv1/vm
cp1=vv1/vp
cm3=0.d0
cp3=0.d0
endif
sm1=sm1+bm(j)*cm1*e1
sp1=sp1+bm(j)*cp1*e1
sm3=sm3+bm(j)*cm3*e1
110 sp3=sp3+bm(j)*cp3*e1
do 120 j=1,nq
br=sbr(j)
e1=e41(j)
vp=1.d0+br*v
vm=1.d0-br*v
vza=v*v+zal(j)
bb2=1.d0-br*br
if (ll) then
cm2=vza/vm
cp2=vza/vp
cm4=val(j)/vm
cp4=val(j)/vp
else
cm2=1.d0+br*v
cp2=1.d0-br*v
cm4=bb2/vm
cp4=bb2/vp
endif
sm2=sm2+bm(nq+j)*cm2*e1
sp2=sp2+bm(nq+j)*cp2*e1
sm4=sm4+bm(nq+j)*cm4*e1
sp4=sp4+bm(nq+j)*cp4*e1
120 continue
a1=bm(1)
aiul(i)=a1+tau+v+sp1+sp2
aidl(i)=a1+tau-v+sm1+sm2
aiur(i)=a1+tau+v+sp3+sp4
aidr(i)=a1+tau-v+sm3+sm4
100 continue
return
end
subroutine mint(aiul,aidl,aiur,aidr,tau,us,ni,zal,val,zbe,vbe,
& da0,da2,da4,db0,db2,db4,bm,qs1,sbr,sbl,nq)
implicit real*8(a-h,o-z)
parameter (nqmax=200)
c this code finds the radiation field for a 2*2 Rayleigh-Cabannes
c scattering homogeneous plane-parallel atmosphere at the angle
c points us(i), (i=1,...,ni).
c
c The conservative Milne problem.
c
c Viik 89 12 10
logical ll
dimension aiul(ni),aidl(ni),aiur(ni),aidr(ni),zal(1),val(1),
& zbe(1),vbe(1),da0(1),da2(1),da4(1),db0(1),db2(1),
& db4(1),el(nqmax),er(nqmax),us(ni),bm(1),
& sbr(nq),sbl(nq)
300 format(1x,i4,2x,3e16.8)
301 format(1x,65('-'))
ll=qs1.ge.(1.d0-1.d-12)
aks=.375d0*qs1
dze=.25d0*(1.d0-qs1)
do 140 i=1,nq
el(i)=dexp(-sbl(i)*tau)
140 er(i)=dexp(-sbr(i)*tau)
c print 301
c do 1 j=1,ni
c print 300,j,us(j)
c 1 continue
c print 301
do 100 j=1,ni
v=us(j)
v2=v*v
vm=1.d0-v2
c print 300,j,v
c if (v.ne.0d0) then
ev=dexp(-tau/v)
c else if (v.eq.0d0) then
c ev=0d0
c endif
x1=bm(1)+v+tau
x2=bm(1)*(1.d0-ev)+tau-v+v*ev
s1=0.d0
s2=0.d0
s3=0.d0
s4=0.d0
s5=0.d0
s6=0.d0
s7=0.d0
s8=0.d0
s9=0.d0
s10=0.d0
s11=0.d0
s12=0.d0
do 110 k=1,nq
bmk=bm(k)
cmk=bm(k+nq)
b4=sbr(k)
b5=sbl(k)
upb=1.d0+b5*v
umb=1.d0-b5*v
upa=1.d0+b4*v
uma=1.d0-b4*v
eb=el(k)
ea=er(k)
if (ll) then
s1=s1+bmk*eb/upb
s2=s2+cmk*uma*ea
s3=s3+bmk*(eb-ev)/umb
s4=s4+cmk*upa*(ea-ev)
s6=s6+cmk*(1.d0-b4*b4)*ea/upa
s12=s12+cmk*(1.d0-b4*b4)*(ea-ev)/uma
else
pa1=da2(k)-da4(k)+zal(k)*(da0(k)-da2(k))
pb1=db2(k)-db4(k)+zbe(k)*(db0(k)-db2(k))
pa2=da2(k)+da0(k)*(zal(k)+val(k))
pb2=db2(k)+db0(k)*(zbe(k)+vbe(k))
ga1=da4(k)+zal(k)*da2(k)+val(k)*da0(k)
gb1=db4(k)+zbe(k)*db2(k)+vbe(k)*db0(k)
s1=s1+bmk*pb1*eb/upb
s2=s2+cmk*pa1*ea/upa
s3=s3+bmk*pb2*eb/upb
s4=s4+cmk*pa2*ea/upa
s5=s5+bmk*gb1*eb/upb
s6=s6+cmk*ga1*ea/upa
s7=s7+bmk*pb1*(eb-ev)/umb
s8=s8+cmk*pa1*(ea-ev)/uma
s9=s9+bmk*pb2*(eb-ev)/umb
s10=s10+cmk*pa2*(ea-ev)/uma
s11=s11+bmk*gb1*(eb-ev)/umb
s12=s12+cmk*ga1*(ea-ev)/uma
endif
if (k.eq.1) then
s1=0.d0
s3=0.d0
s5=0.d0
s7=0.d0
s9=0.d0
s11=0.d0
endif
110 continue
if (.not.ll) then
alp=x1+2.d0*aks*(s1+s2)+dze*(s3+s4)
arp=x1+aks*(s5+s6)+dze*(s3+s4)
alm=x2+2.d0*aks*(s7+s8)+dze*(s9+s10)
arm=x2+aks*(s11+s12)+dze*(s9+s10)
aiul(j)=(1.d0-v2)*alp+v2*arp
aidl(j)=(1.d0-v2)*alm+v2*arm
aiur(j)=arp
aidr(j)=arm
else
aiul(j)=x1+vm*s1+s2
aidl(j)=x2+vm*s3+s4
aiur(j)=x1+s6
aidr(j)=x2+s12
endif
100 continue
return
end
subroutine mssqq(si,sq,qi,qq,ql,qr,tau,qs1,sbr,sbl,zal,val,
& zbe,vbe,da0,da2,da4,db0,db2,db4,bm,nq)
implicit real*8(a-h,o-z)
c this code finds the source functions s_i and s_q together with
c q_i and q_q functions for the 2*2 Rayleigh-Cabannes
c scattering homogeneous plane-parallel atmosphere
c
c The conservative Milne problem
c
c Viik 92 05 12
dimension zal(nq),val(nq),zbe(nq),vbe(nq),da0(nq),da2(nq),
& da4(nq),db0(nq),db2(nq),db4(nq),bm(1),
& sbr(nq),sbl(nq)
qs2=1d0-qs1
if (qs1.ne.1d0) then
s2=0.d0
s4=0.d0
s6=0.d0
do 10 k=1,nq
cmk=bm(k+nq)
er=dexp(-sbr(k)*tau)
pa2=da2(k)+da0(k)*(zal(k)+val(k))
pa1=da2(k)-da4(k)+zal(k)*(da0(k)-da2(k))
ga1=pa2-pa1
s2=s2+cmk*pa1*er
s4=s4+cmk*pa2*er
s6=s6+cmk*ga1*er
10 continue
s1=0.d0
s3=0.d0
s5=0.d0
do 20 k=2,nq
bmk=bm(k)
el=dexp(-sbl(k)*tau)
pb2=db2(k)+db0(k)*(zbe(k)+vbe(k))
pb1=db2(k)-db4(k)+zbe(k)*(db0(k)-db2(k))
gb1=pb2-pb1
s1=s1+bmk*pb1*el
s3=s3+bmk*pb2*el
s5=s5+bmk*gb1*el
20 continue
ql=bm(1)+0.75d0*qs1*(s1+s2)+0.25d0*qs2*(s3+s4)
qr=bm(1)+0.375d0*qs1*(s5+s6)+0.25d0*qs2*(s3+s4)
else
s2=0d0
s3=0d0
do 30 k=1,nq
za=sbr(k)
dex=dexp(-za*tau)
s2=s2+bm(k+nq)*dex
s3=s3+bm(k+nq)*(1d0-za*za)*dex
30 continue
s1=0d0
do 40 k=2,nq
s1=s1+bm(k)*dexp(-sbl(k)*tau)
40 continue
ql=bm(1)+s1+s2
qr=bm(1)+s3
endif
qi=2d0/3d0*(ql+2d0*qr)
qq=2d0/3d0*dsqrt(2d0/qs1)*(ql-qr)
sl=tau+ql
sr=tau+qr
si=2d0/3d0*(sl+2d0*sr)
sq=2d0/3d0*dsqrt(2d0/qs1)*(sl-sr)
return
end
c
c
c
c NONCONSERVATIVE CASE
c
c
c
subroutine wurzel(z,vig,alamb,alamq)
implicit real*8(a-h,o-z)
c
c this code finds the solution of the characteristic
c equation for vector transfer
c
c Viik 92 02 28
external rfx
tol=1d-14
qol=1d0-tol
if (alamb.ge.qol.or.alamq.ge.qol) then
z=0d0
vig=0d0
return
endif
if (alamb.le.tol.and.alamq.le.tol) then
z=1d0
vig=0d0
return
endif
f7=alamb+alamq-.1d0
alq=alamb*alamq
if (f7.le.0d0) then
d7=(28d0-45d0*alq)/(14d0*alamb+1d1*alamq-3d1*alq)
z=1d0-2d0*dexp(-d7)
vig=0d0
return
endif
ax=0d0
bx=1d0-tol
npun=500
x0=ax
jmax=65
y0=rfx(ax,alamb,alamq)
do 20 i=1,npun
x=bx/dfloat(npun)*dfloat(i)
y=rfx(x,alamb,alamq)
if (y*y0.lt.0d0) then
z=rtbisx(x0,x,rfx,tol,noz,jmax,alamb,alamq)
vig=rfx(z,alamb,alamq)
return
endif
y0=y
x0=x
c print 306,alamb,alamq,x,y
20 continue
return
end
function rfx(x,alamb,alamq)
implicit real*8(a-h,o-z)
c this code finds the left hand side of the characteristic
c equation for the Rayleigh-Cabannes scattering
c Viik 90 01 23
f0=(1d0-alamb)*(1d0-alamq)
if (x.eq.0d0) then
rfx=112d0*f0
else
x2=x*x
x3=x2*x
x4=x3*x
if (alamb.lt.0.999999d0) then
g0=dlog((1d0+x)/(1d0-x))/x
g2=1d0/x2*(-2d0+g0)
g4=-2d0/x2*(1d0/x2+1d0/3d0)+1d0/x4*g0
else
g0=2d0+2d0/3d0*x2+2d0/5d0*x4+2d0/7d0*x2*x4+2d0/9d0*
& x4*x4+2d0/11d0*x2*x4*x4
g2=2d0/3d0+2d0/5d0*x2+2d0/7d0*x4+2d0/9d0*x2*x4+
& 2d0/11d0*x4*x4+2d0/13d0*x2*x4*x4
g4=2d0/5d0+2d0/7d0*x2+2d0/9d0*x4+2d0/11d0*x2*x4+
& 2d0/13d0*x4*x4+2d0/15d0*x2*x4*x4
endif
y1=alamb*(g0-g2)
y2=alamq*(g0-g2)
rfx=(3d0*y1-4d0)*(15d0*y2-28d0)+
& 4d0*(7d0*alamb-1d1*alamq)*g0+
& 12d0*(5d0*alamq*(3d0-alamb)-7d0*alamb)*g2-
& 180d0*alamq*(1d0-alamb)*g4
endif
return
end
function rtbisx(x1,x2,func,xacc,noz,jmax,alamb,alamq)
implicit real*8(a-h,o-z)
c using bisection, find the root of a function func known
c to lie between x1 and x2. the root, returned as rtbis
c will be refined jmax times until its accuracy is xacc.
c alamb and alamq are parameters.
300 format(1x,'no zeros or even number of zeros')
noz=0
fmid=func(x2,alamb,alamq)
f=func(x1,alamb,alamq)
if (f*fmid.gt.0.d0) then
noz=1
print 300
return
endif
if (f.lt.0.d0) then
rtbisx=x1
dx=x2-x1
else
rtbisx=x2
dx=x1-x2
endif
do 11 j=1,jmax
dx=dx*.5d0
xmid=rtbisx+dx
fmid=func(xmid,alamb,alamq)
if (fmid.lt.0.d0) rtbisx=xmid
if (dabs(dx).lt.xacc.or.fmid.eq.0.d0) return
11 continue
return
end
subroutine zava(zk,vk,zal,val,zbe,vbe,alamb,x,
& qs1,sbr,sbl,uu,cc,nq)
implicit real*8(a-h,o-z)
c this code finds for Rayleigh scattering
c the values of functions z(k) and v(k)
c the nonconservative Milne scattering
c Viik 92 05 03
dimension zal(1),val(1),zbe(1),vbe(1),sbr(nq),sbl(nq),
& uu(nq),cc(nq)
qs2=1.d0-qs1
aks=.375d0*alamb*qs1
dze=.25d0*alamb*qs2
do 10 ix=1,2
do 30 i=1,nq
be=sbr(i)
if (ix.eq.2) be=sbl(i)
d0=0.d0
d2=0.d0
d4=0.d0
do 40 j=1,nq
u2=uu(j)*uu(j)
u4=u2*u2
be2=u2*be*be
cj=cc(j)/(1.d0-be2)
d0=d0+cj
d2=d2+cj*u2
d4=d4+cj*u4
40 continue
d0=2d0*d0
d2=2d0*d2
d4=2d0*d4
a11=1.d0-2.d0*aks*(d0-d2)-dze*d0
a12=-dze*d0
a21=-aks*d2-dze*d0
a22=1.d0-d0*(aks+dze)
b1=2.d0*aks*(d2-d4)+dze*d2
b2=aks*d4+dze*d2
del=a11*a22-a12*a21
z=(b1*a22-b2*a12)/del
v=(b2*a11-b1*a21)/del
if (ix.ne.1) go to 50
zal(i)=z
val(i)=v
go to 30
50 zbe(i)=z
vbe(i)=v
30 continue
10 continue
dd0=0.d0
dd2=0.d0
dd4=0.d0
do 60 j=1,nq
u2=uu(j)*uu(j)
u4=u2*u2
be2=u2*x*x
cj=cc(j)/(1.d0-be2)
dd0=dd0+cj
dd2=dd2+cj*u2
dd4=dd4+cj*u4
60 continue
dd0=2d0*dd0
dd2=2d0*dd2
dd4=2d0*dd4
a11=1d0-2d0*aks*(dd0-dd2)-dze*dd0
a12=-dze*dd0
a21=-aks*dd2-dze*dd0
a22=1d0-dd0*(aks+dze)
b1=2d0*aks*(dd2-dd4)+dze*dd2
b2=aks*dd4+dze*dd2
del=a11*a22-a12*a21
zk=(b1*a22-b2*a12)/del
vk=(b2*a11-b1*a21)/del
return
end
subroutine coefmi(bm,zal,val,zbe,vbe,zk,vk,x,am,sbr,sbl,uu,nq)
implicit real*8(a-h,o-z)
c this code finds the coefficients a(i) and c(i)
c in a semi-infinite homogeneous atmosphere
c for 2*2 Rayleigh scattering
c
c the Milne problem
c
c the scattering is not conservative!
c
c
c Viik 92 05 03
c
300 format(1x,8e15.5)
301 format(1x,120(1h-))
dimension zal(nq),val(nq),zbe(nq),vbe(nq),
& am(2*nq*2*nq),bm(2*nq),sbr(nq),sbl(nq),uu(nq)
n2=2*nq
do 10 j=1,n2
do 20 i=1,n2
ij=i+n2*(j-1)
if (i.le.nq.and.j.le.nq) then
am(ij)=(uu(i)*uu(i)+zbe(j))/(1d0-uu(i)*sbl(j))
else if (i.le.nq.and.j.gt.nq) then
am(ij)=(uu(i)*uu(i)+zal(j-nq))/
& (1d0-uu(i)*sbr(j-nq))
else if (i.gt.nq.and.j.le.nq) then
am(ij)=vbe(j)/(1d0-uu(i-nq)*sbl(j))
else if (i.gt.nq.and.j.gt.nq) then
am(ij)=val(j-nq)/(1d0-uu(i-nq)*sbr(j-nq))
endif
20 continue
if (j.le.nq) then
bm(j)=-(zk+uu(j)*uu(j))/(1d0+x*uu(j))
else
bm(j)=-vk/(1d0+x*uu(j-nq))
endif
10 continue
call linsys(am,bm,n2,mis)
return
end
subroutine mncint(aiul,aidl,aiur,aidr,alamb,tau,us,ni,
& x,qs1,zk,vk,zal,val,zbe,vbe,
& da0,da2,da4,db0,db2,db4,dk0,dk2,dk4,bm,
& sbr,sbl,nq)
implicit real*8(a-h,o-z)
c this code finds the radiation field for a 2*2 Rayleigh-Cabannes
c scattering homogeneous plane-parallel atmosphere at the angle
c points us(i), (i=1,...,ni).
c the nonconservative Milne problem.
c Viik 92 05 12
dimension aiul(1),aidl(1),aiur(1),aidr(1),zal(1),val(1),
& zbe(1),vbe(1),da0(1),da2(1),da4(1),db0(1),db2(1),
& db4(1),el(48),er(48),us(1),bm(1),sbr(nq),sbl(nq)
aks=.375d0*alamb*qs1
dze=.25d0*alamb*(1.d0-qs1)
dex=dexp(x*tau)
do 140 i=1,nq
el(i)=dexp(-sbl(i)*tau)
140 er(i)=dexp(-sbr(i)*tau)
do 100 j=1,ni
v=us(j)
v2=v*v
vm=1.d0-v2
ev=dexp(-tau/v)
gkm=dex/(1d0-x*v)
gkp=(dex-ev)/(1d0+x*v)
pk1=dk2-dk4+zk*(dk0-dk2)
pk2=dk2+dk0*(zk+vk)
gk1=dk4+zk*dk2+vk*dk0
s1=0.d0
s2=0.d0
s3=0.d0
s4=0.d0
s5=0.d0
s6=0.d0
s7=0.d0
s8=0.d0
s9=0.d0
s10=0.d0
s11=0.d0
s12=0.d0
do 110 k=1,nq
bmk=bm(k)
cmk=bm(k+nq)
b4=sbr(k)
b5=sbl(k)
upb=1.d0+b5*v
umb=1.d0-b5*v
upa=1.d0+b4*v
uma=1.d0-b4*v
pa1=da2(k)-da4(k)+zal(k)*(da0(k)-da2(k))
pb1=db2(k)-db4(k)+zbe(k)*(db0(k)-db2(k))
pa2=da2(k)+da0(k)*(zal(k)+val(k))
pb2=db2(k)+db0(k)*(zbe(k)+vbe(k))
ga1=da4(k)+zal(k)*da2(k)+val(k)*da0(k)
gb1=db4(k)+zbe(k)*db2(k)+vbe(k)*db0(k)
eb=el(k)
ea=er(k)
s1=s1+bmk*pb1*eb/upb
s2=s2+cmk*pa1*ea/upa
s3=s3+bmk*pb2*eb/upb
s4=s4+cmk*pa2*ea/upa
s5=s5+bmk*gb1*eb/upb
s6=s6+cmk*ga1*ea/upa
s7=s7+bmk*pb1*(eb-ev)/umb
s8=s8+cmk*pa1*(ea-ev)/uma
s9=s9+bmk*pb2*(eb-ev)/umb
s10=s10+cmk*pa2*(ea-ev)/uma
s11=s11+bmk*gb1*(eb-ev)/umb
s12=s12+cmk*ga1*(ea-ev)/uma
110 continue
alp=2d0*aks*(s1+s2+pk1*gkm)+dze*(s3+s4+pk2*gkm)
arp=aks*(s5+s6+gk1*gkm)+dze*(s3+s4+pk2*gkm)
alm=2d0*aks*(s7+s8+pk1*gkp)+dze*(s9+s10+pk2*gkp)
arm=aks*(s11+s12+gk1*gkp)+dze*(s9+s10+pk2*gkp)
aiul(j)=(1.d0-v2)*alp+v2*arp
aidl(j)=(1.d0-v2)*alm+v2*arm
aiur(j)=arp
aidr(j)=arm
100 continue
return
end
subroutine mintgn(aul,adl,aur,adr,tau,x,bm,zk,vk,zal,val,
& zbe,vbe,uu,sbr,sbl,nq)
implicit real*8(a-h,o-z)
c
c this code finds intensities at any level for
c 2*2 Rayleigh-Cabannes scattering in homogeneous
c plane-parallel non-conservative atmosphere at
c Gaussian points uu(i), i=1,...,nq.
c
c the Milne problem
c
c Viik 02 10 98
c
dimension bm(1),aul(nq),adl(nq),
& aur(nq),adr(nq),zal(nq),val(nq),zbe(nq),vbe(nq),
& sbr(nq),sbl(nq),uu(nq)
d=dexp(x*tau)
do 20 i=1,nq
u=uu(i)
s1=0d0
s2=0d0
s3=0d0
s4=0d0
do 10 k=1,nq
e5=dexp(-tau*sbl(k))
e4=dexp(-tau*sbr(k))
qpbe=(u*u+zbe(k))/(1d0+u*sbl(k))
qmbe=(u*u+zbe(k))/(1d0-u*sbl(k))
qpal=(u*u+zal(k))/(1d0+u*sbr(k))
qmal=(u*u+zal(k))/(1d0-u*sbr(k))
rpbe=vbe(k)/(1d0+u*sbl(k))
rmbe=vbe(k)/(1d0-u*sbl(k))
rpal=val(k)/(1d0+u*sbr(k))
rmal=val(k)/(1d0-u*sbr(k))
s1=s1+bm(k)*qpbe*e5+bm(k+nq)*qpal*e4
s2=s2+bm(k)*rpbe*e5+bm(k+nq)*rpal*e4
s3=s3+bm(k)*qmbe*e5+bm(k+nq)*qmal*e4
s4=s4+bm(k)*rmbe*e5+bm(k+nq)*rmal*e4
10 continue
ausl=(zk+u*u)/(1d0-x*u)*d
ausr=vk/(1d0-x*u)*d
aul(i)=s1+ausl
aur(i)=s2+ausr
adl(i)=s3+(zk+u*u)/(1d0+x*u)*d
adr(i)=s4+vk/(1d0+x*u)*d
20 continue
return
end
subroutine mncqq(qi,qq,alamb,tau,qs1,zk,vk,sbr,sbl,
& zal,val,zbe,vbe,da0,da2,da4,db0,db2,db4,
& dk0,dk2,dk4,bm,nq)
implicit real*8(a-h,o-z)
c this code finds the q_i and q_q functions
c for the 2*2 Rayleigh-Cabannes scattering homogeneous
c plane-parallel atmosphere
c
c the nonconservative Milne problem
c
c Viik 92 05 12
dimension zal(1),val(1),zbe(1),vbe(1),da0(1),da2(1),
& da4(1),db0(1),db2(1),db4(1),bm(1),
& sbr(nq),sbl(nq)
aks=.375d0*alamb*qs1
dze=.25d0*alamb*(1.d0-qs1)
s1=0.d0
s2=0.d0
s3=0.d0
s4=0.d0
s5=0.d0
s6=0.d0
pk1=dk2-dk4+zk*(dk0-dk2)
pk2=dk2+(zk+vk)*dk0
gk1=dk4+zk*dk2+vk*dk0
do 10 k=1,nq
bmk=bm(k)
cmk=bm(k+nq)
el=dexp(-sbl(k)*tau)
er=dexp(-sbr(k)*tau)
pb2=db2(k)+db0(k)*(zbe(k)+vbe(k))
pa2=da2(k)+da0(k)*(zal(k)+val(k))
pb1=db2(k)-db4(k)+zbe(k)*(db0(k)-db2(k))
pa1=da2(k)-da4(k)+zal(k)*(da0(k)-da2(k))
gb1=db4(k)+zbe(k)*db2(k)+vbe(k)*db0(k)
ga1=da4(k)+zal(k)*da2(k)+val(k)*da0(k)
s1=s1+bmk*pb1*el
s2=s2+cmk*pa1*er
s3=s3+bmk*pb2*el
s4=s4+cmk*pa2*er
s5=s5+bmk*gb1*el
s6=s6+cmk*ga1*er
10 continue
dq=dze*(s3+s4)
ql=2d0*aks*(s1+s2)+dq
qr=aks*(s5+s6)+dq
qi=2d0/3d0*(ql+2d0*qr)
qq=2d0/3d0*dsqrt(2d0/qs1)*(ql-qr)
return
end
subroutine mncsou(soul,sour,si,sq,qi,qq,alamb,tau,qs1,x,zk,vk,
& zal,val,zbe,vbe,da0,da2,da4,db0,db2,db4,
& dk0,dk2,dk4,sbr,sbl,bm,nq)
implicit real*8(a-h,o-z)
c this code finds the s sub l and s sub r functions together
c with s_i and s_q
c for the 2*2 Rayleigh-Cabannes scattering homogeneous
c plane-parallel atmosphere
c the nonconservative Milne problem
c Viik 92 05 12
dimension zal(nq),val(nq),zbe(nq),vbe(nq),da0(nq),da2(nq),
& da4(nq),db0(nq),db2(nq),db4(nq),bm(2*nq),
& sbr(nq),sbl(nq)
aks=.375d0*alamb*qs1
dze=.25d0*alamb*(1.d0-qs1)
s1=0.d0
s2=0.d0
s3=0.d0
s4=0.d0
s5=0.d0
s6=0.d0
dex=dexp(tau*x)
pk1=dk2-dk4+zk*(dk0-dk2)
pk2=dk2+(zk+vk)*dk0
gk1=dk4+zk*dk2+vk*dk0
do 10 k=1,nq
bmk=bm(k)
cmk=bm(k+nq)
el=dexp(-sbl(k)*tau)
er=dexp(-sbr(k)*tau)
pb2=db2(k)+db0(k)*(zbe(k)+vbe(k))
pa2=da2(k)+da0(k)*(zal(k)+val(k))
pb1=db2(k)-db4(k)+zbe(k)*(db0(k)-db2(k))
pa1=da2(k)-da4(k)+zal(k)*(da0(k)-da2(k))
gb1=db4(k)+zbe(k)*db2(k)+vbe(k)*db0(k)
ga1=da4(k)+zal(k)*da2(k)+val(k)*da0(k)
s1=s1+bmk*pb1*el
s2=s2+cmk*pa1*er
s3=s3+bmk*pb2*el
s4=s4+cmk*pa2*er
s5=s5+bmk*gb1*el
s6=s6+cmk*ga1*er
10 continue
dq=dze*(s3+s4)
de=dze*(s3+s4+pk2*dex)
soul=2d0*aks*(s1+s2+pk1*dex)+de
sour=aks*(s5+s6+gk1*dex)+de
ql=2d0*aks*(s1+s2)+dq
qr=aks*(s5+s6)+dq
si=2d0/3d0*(soul+2d0*sour)
sq=2d0/3d0*dsqrt(2d0/qs1)*(soul-sour)
qi=2d0/3d0*(ql+2d0*qr)
qq=2d0/3d0*dsqrt(2d0/qs1)*(ql-qr)
return
end