c
c A PACKAGE OF SUBROUTINES FOR SOLVING THE VECTOR EQUATION OF
c TRANSFER IN A PLANETARY ATMOSPHERE WHERE THE SCATTERING OBEYS
c THE RAYLEIGH-CABANNES LAW.
c
c The theoretical background is given in:
c
c 1. Viik T., Rayleigh Scattering in a Homegeneous Planeparallel
c Atmosphere, Tallinn, Valgus, 1989 (alas, in Russian)
c
c 2. Viik T., Rayleigh Scattering in Planetary Atmospheres,
c Earth, Moon, and Planets, vol. 46, pp.261-284, 1989.
c
c The main input parameters are:
c
c tau0 - optical thickness of the atmosphere
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 uinc - cosine of the angle between the direction of incident
c beam and the normal to the atmosphere (uinc=1 for
c normal incidence and uinc=0 for grazing angle)
c us - array of cosines of arbitrary angular points
c us(1,...,ni)
c ni - dimension of array us
c alb - albedo of the Lambert bottom of the atmosphere
c nq - order of the Gaussian quadrature
c uu - array of the Gaussian quadrature points
c uu(1,...,nq)
c cc - array of the Gaussian quadrature weights
c cc(1,...,nq)
c FL, FR, FU, FV - the components of the Stokes vector of the
c incident beam
c tau - optical thickness, measured from the upper boundary
c of the atmosphere
c fi - azimuth angle
c
c For checking purposes I have included the code
c CHECK_UPA.FOR
c together with a file of results
c RESULT_RAYLEIGH.DAT.
c
c
c
c GENERAL DESCRIPTION
c
c GAULEG - finds the array of points x(1...nq) and the array of weights
c w(1...nq) 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 ALBEGA - finds for Rayleigh-Cabannes scattering the values of
c functions alpha, beta and gamma
c
c COEFCON - finds the coefficients a(i), b(i), c(i) and d(i)
c for the 2*2 Rayleigh-Cabannes scattering in a
c conservative atmosphere
c
c COEFNCON - finds the coefficients a(i), b(i), c(i) and d(i)
c for the 2*2 Rayleigh-Cabannes scattering in a
c non-conservative atmosphere
c
c INTGAUSSCON - finds the intensities at arbitrary optical depth
c for the 2*2 Rayleigh-Cabannes scattering
c at Gaussian points in conservative atmosphere
c (l,r) representation
c
c INTGAUSSNCON - finds the intensities at arbitrary optical depth
c for the 2*2 Rayleigh-Cabannes scattering
c at Gaussian points in non-conservative atmosphere
c (l,r) representation
c
c INTARBLR - finds the intensities at arbitrary optical depth
c for the 2*2 Rayleigh-Cabannes scattering
c at arbitrary angular points us(1,...,ni)
c (l,r) representation
c
c GAUSSINT - finds the Stokes vector of intensity at arbitrary
c optical depth for the 4*4 Rayleigh-Cabannes scattering
c at Gaussian angular points
c Both (l,r) and (I,Q) representations
c
c ALLINT - finds the Stokes vector of intensity at arbitrary
c optical depth for the 4*4 Rayleigh-Cabannes scattering
c at arbitrary angular points us(1,...,ni)
c Both (l,r) and (I,Q) representations
c
c STOKES - finds the Stokes vector of intensity at arbitrary
c optical depth for the 4*4 Rayleigh scattering
c at arbitrary angular points us(1,...,ni).
c ONLY THE CONSERVATIVE CASE OF THE RAYLEIGH
c SCATTERING IS CONSIDERED BY CHANDRASEKHAR METHOD.
c BY USING STOKES YOU MAY FIND A MORE ACCURATE VERSION
c OF COULSON-DAVE-SEKERA TABLES!
c
c
c CODES FOR ISOTROPICALLY SCATTERING ATMOSPHERE
c NEEDED IN VECTOR CASE
c
c SICOEFF - finds the coefficients of the Sobolev resolvent
c function for optically semi-infinite atmosphere
c
c HFUN - finds the H-function
c
c HIGIF - finds the h- and g-functions
c
c DGKAP - finds the derivatives of h- and g-functions with
c respect to angular variable
c
c SINTST - finds the intensities in an optically semi-infinite
c atmosphere at arbitrary angular points us(1,...,ni)
c
c FICOEFF - finds the coefficients of the Sobolev resolvent
c function for optically infinite atmosphere
c (planetary problem)
c
c XYFUN - finds the X- and Y-functions
c
c XIYIF - finds the x- and y-functions
c
c DXIDYI - finds the derivatives of x- and y-functions with
c respect to angular variable
c
c FINTST - finds the intensities in an optically infinite
c atmosphere at arbitrary angular points us(1,...,ni)
c (planetary problem)
c
c
c GENERAL CODES
c
c
c
subroutine gauleg(x1,x2,x,w,nq)
implicit real*8(a-h,o-z)
c
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 nq
c containing the abscissas and weigths of the gauss-
c legendre nq-point quadrature formula.
c A code by G.B.Rybicki from NUMERICAL RECIPES by W.H.Press,
c B.P.Flannery, S.A.Teukolsky, and W.T.Vetterling
c
parameter (eps=.11d-15,pi=3.141592653589793238d0)
dimension x(nq),w(nq)
m=(nq+1)/2
itmax=40
xm=.5d0*(x2+x1)
xl=.5d0*(x2-x1)
do 10 i=1,m
z =dcos(pi*(dfloat(i)-.25d0)/(dfloat(nq)+.5d0))
it=0
30 p1=1.d0
p2=0.d0
do 20 j=1,nq
dj=dfloat(j)
p3=p2
p2=p1
p1=((2.d0*dj-1.d0)*z*p2-(dj-1.d0)*p3)/dj
20 continue
pp=dfloat(nq)*(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 30
x(i)=xm-xl*z
x(nq+1-i)=xm+xl*z
w(i)=2.d0*xl/((1.d0-z*z)*pp*pp)
w(nq+1-i)=w(i)
10 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
c
c linsys - subroutine for solving a system
c of linear algebraic equations,
c e.g. - system/360 scientific subroutine package
c (360a-cm-03x) version 3.
c 4th edition, IBM, Technical Publications Dept.,
c N.Y., 1960-1970
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=0
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
function charfunct(iw,u,alamb)
implicit real*8(a-h,o-z)
c charfunct is the characteristic function for
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 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
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
function apprchaf(x,alamb,qs1,uu,cc,nq)
implicit real*8(a-h,o-z)
c vt858 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
subroutine zeropol(sbr,sbl,alamb,qs1,uu,cc,nq)
implicit real*8(a-h,o-z)
parameter (nmax=400,siu=1d-10,nqmax=200)
c this code evaluates the zeros of the approximate
c characteristic equation for Rayleigh 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=0.d0
bx=u-siu
if (ij.ne.1) ax=1.d0/uu(ji+1)+siu
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 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 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 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,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
subroutine albega(zal,val,zbe,vbe,alamb,uinc,qs1,fl,fr,
& aa,bb,kp2,kp,sbr,sbl,uu,cc,nq)
implicit real*8(a-h,o-z)
c this code finds for Rayleigh scattering
c the values of
c alfa,beta and gamma (in array bb))
c and functions z(k) and v(k)
c Viik 87 12 22
logical logic
dimension aa(kp2),bb(kp),zal(1),val(1),zbe(1),vbe(1),
& sbr(nq),sbl(nq),uu(nq),cc(nq)
stol=1.d0-1.d-14
qs2=1.d0-qs1
aks=.375d0*alamb*qs1
dze=.25d0*alamb*qs2
logic=alamb.ge.stol.and.qs1.ge.stol
if (.not.logic) then
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
be3=1.d0-be2
cj=cc(j)/be3
d0=d0+cj
d2=d2+cj*u2
40 d4=d4+cj*u4
d0=2.d0*d0
d2=2.d0*d2
d4=2.d0*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.eq.1) then
zal(i)=z
val(i)=v
else
zbe(i)=z
vbe(i)=v
endif
30 continue
10 continue
endif
d0=0.d0
do 60 i=1,nq
u=uu(i)/uinc
uh=1.d0-u*u
60 d0=d0+cc(i)/uh
d0=2.d0*d0
uin2=uinc*uinc
d2=uin2*(d0-2.d0)
d4=uin2*uin2*(d0-2.d0)-2.d0*uin2/3.d0
aa(1)=1.d0-aks*(3.d0*d4-2.d0*d2)
aa(2)=-2.d0*aks*(d2-d4)-dze*d2
aa(3)=-aks*d4-dze*d2
aa(4)=-aks*(3.d0*d2-2.d0*d0)
aa(5)=1.d0-2.d0*aks*(d0-d2)-dze*d0
aa(6)=-aks*d2-dze*d0
aa(7)=-aks*d0
aa(8)=-dze*d0
aa(9)=1.d0-d0*(aks+dze)
bb(1)=.5d0*aks*((3.d0*uin2-2.d0)*fl+fr)
bb(2)=aks*(1.d0-uin2)*fl+.5d0*dze*(fl+fr)
bb(3)=.5d0*aks*(uin2*fl+fr)+.5d0*dze*(fl+fr)
call linsys(aa,bb,3,ks)
return
end
subroutine coefcon(bm,tau0,uinc,alb,fl,fr,bb,kp,am,
& sbr,sbl,uu,cc,nq)
implicit real*8(a-h,o-z)
parameter (nqmax=200)
c this code finds for Rayleigh (2*2) problem
c the coefficients a(i),b(i),c(i) and d(i)
c if (alamb.eq.1.and.qs1.eq.1)
c Viik 87 12 22
c
c
dimension jmin(16),jmax(16),imin(16),imax(16),
& bb(kp),am(1),bm(1),exl(nqmax),exr(nqmax),
& dpb(nqmax),dmb(nqmax),dpa(nqmax),dma(nqmax),
& sbr(nq),sbl(nq),uu(nq),cc(nq)
do 65 i=1,nq
exl(i)=dexp(-tau0*sbl(i))
65 exr(i)=dexp(-tau0*sbr(i))
if (uinc.eq.0.d0) pause "uinc is zero"
ext=dexp(-tau0/uinc)
n4=4*nq
n3=3*nq
n2=2*nq
c if (tau0.gt.1.d5) go to 130
nq1=nq+1
n21=n2+1
n31=n3+1
ij=0
do 99 i=1,4
do 99 j=1,4
ij=ij+1
imin(ij)=(i-1)*nq+1
jmin(ij)=(j-1)*nq+1
imax(ij)=i*nq
99 jmax(ij)=j*nq
if (alb.gt.1.d-12) go to 11
do 100 ki=1,16
j1=jmin(ki)
j2=jmax(ki)
i1=imin(ki)
i2=imax(ki)
do 101 i=i1,i2
ie=i-i1+1
u=uu(ie)
u1=1.d0-u*u
do 101 j=j1,j2
je=j-j1+1
s=sbl(je)
if (j.gt.n2) s=sbr(je)
s1=1.d0-s*s
e=exl(je)
if (j.gt.n2) e=exr(je)
sp=1.d0+u*s
sm=1.d0-u*s
if (sp.eq.0.d0.or.sm.eq.0.d0) pause 'this is gibberish'
ij=i+n4*(j-1)
if (ki.ne.1) go to 102
am(ij)=u1*e/sp
if (je.eq.1) am(ij)=1.d0
102 if (ki.ne.2) go to 103
am(ij)=u1/sm
if (je.eq.1) am(ij)=tau0+u
103 if (ki.eq.3) am(ij)=sm*e
if (ki.eq.4) am(ij)=sp
if (ki.ne.5) go to 104
am(ij)=u1/sm
if (je.eq.1) am(ij)=1.d0
104 if (ki.ne.6) go to 105
am(ij)=u1*e/sp
if (je.eq.1) am(ij)=-u
105 if (ki.eq.7) am(ij)=sp
if (ki.eq.8) am(ij)=sm*e
if (ki.ne.9) go to 106
am(ij)=0.d0
if (je.eq.1) am(ij)=1.d0
106 if (ki.ne.10) go to 107
am(ij)=0.d0
if (je.eq.1) am(ij)=tau0+u
107 if (ki.eq.11) am(ij)=s1*e/sp
if (ki.eq.12) am(ij)=s1/sm
if (ki.ne.13) go to 108
am(ij)=0.d0
if (je.eq.1) am(ij)=1.d0
108 if (ki.ne.14) go to 109
am(ij)=0.d0
if (je.eq.1) am(ij)=-u
109 if (ki.eq.15) am(ij)=s1/sm
if (ki.eq.16) am(ij)=s1*e/sp
101 continue
100 continue
do 110 i=1,n4
if (i.le.nq) bm(i)=-(bb(1)*uu(i)*uu(i)+bb(2))*
& ext/(1.d0+uu(i)/uinc)
if (i.gt.nq.and.i.le.n2) bm(i)=-(bb(1)*uu(i-nq)*uu(i-nq)+
& bb(2))/(1.d0-uu(i-nq)/uinc)
if (i.gt.n2.and.i.le.n3) bm(i)=-bb(3)*ext/
& (1.d0+uu(i-n2)/uinc)
if (i.gt.n3) bm(i)=-bb(3)/(1.d0-uu(i-n3)/uinc)
110 continue
10 call linsys(am,bm,n4,ks)
return
c 130 nq1=nq+1
c ij=0
c do 131 i=1,2
c do 131 j=1,2
c ij=ij+1
c imin(ij)=(i-1)*nq+1
c imax(ij)=i*nq
c jmin(ij)=(j-1)*nq+1
c 131 jmax(ij)=j*nq
c do 132 ki=1,4
c j1=jmin(ki)
c j2=jmax(ki)
c i1=imin(ki)
c i2=imax(ki)
c do 133 i=i1,i2
c ie=i-i1+1
c u=uu(ie)
c u1=1.d0-u*u
c do 133 j=j1,j2
c je=j-j1+1
c ij=i+n2*(j-1)
c if (ki.ne.1) go to 134
c am(ij)=u1*exl(je)/(1.d0-u*sbl(je))
c if (je.eq.1) am(ij)=1.d0
c 134 if (ki.eq.2) am(ij)=1.d0+sbr(je)*u
c if (ki.ne.3) go to 135
c am(ij)=0.d0
c if (je.eq.1) am(ij)=1.d0
c 135 if (ki.eq.4) am(ij)=(1.d0-sbr(je)*sbr(je))/
c & (1.d0-sbr(je)*u)
c 133 continue
c 132 continue
c do 136 i=1,n2
c if (i.le.nq) bm(i)=-(bb(1)*uu(i)*uu(i)+bb(2))/
c & (1.d0-uu(i)/uinc)
c if (i.gt.nq) bm(i)=-bb(3)/(1.d0-uu(i-nq)/uinc)
c 136 continue
c call linsys(am,bm,n2,kes)
11 pi=3.14159265358979d0
al=alb
uw=1.d0/3.d0
xl=1.d0-al
do 30 j=1,nq
s1=0.d0
s2=0.d0
s3=0.d0
s4=0.d0
z2=sbl(j)
z3=sbr(j)
do 20 i=1,nq
z5=z2*uu(i)
z4=z3*uu(i)
s1=s1+cc(i)/(1.d0+z5)
s2=s2+cc(i)/(1.d0-z5)
s3=s3+cc(i)/(1.d0+z4)
s4=s4+cc(i)/(1.d0-z4)
20 continue
x2=1.d0
if (z2.ge.1.d-12) x2=1.d0/z2
x3=1.d0/z3
dbp1=(1.d0-s1)*x2
dbp3=x2*(uw-x2*(.5d0-dbp1))
dap1=(1.d0-s3)*x3
dap3=x3*(uw-x3*(.5d0-dap1))
dbm1=(s2-1.d0)*x2
dbm3=x2*(x2*(dbm1-.5d0)-uw)
dam1=(s4-1.d0)*x3
dam3=x3*(x3*(dam1-.5d0)-uw)
dpb(j)=al*(dbp3-dbp1)
dmb(j)=al*(dbm3-dbm1)
dpa(j)=al*(-.5d0+uw*z3+(z3*z3-1.d0)*dap1)
dma(j)=al*(-.5d0-uw*z3+(z3*z3-1.d0)*dam1)
30 continue
s5=0.d0
do 21 i=1,nq
21 s5=s5+cc(i)/(1.d0-uu(i)/uinc)
dm1=(s5-1.d0)*uinc
dm3=uinc*(uinc*(dm1-.5d0)-uw)
dmm=(bb(1)*dm3+(bb(2)+bb(3))*dm1)*al
do 40 ki=1,16
j1=jmin(ki)
j2=jmax(ki)
i1=imin(ki)
i2=imax(ki)
do 50 i=i1,i2
ie=i-i1+1
u=uu(ie)
u1=1.d0-u*u
do 50 j=j1,j2
je=j-j1+1
s=sbl(je)
w=dpb(je)
t=dmb(je)
e=exl(je)
if (j.le.n2) go to 60
s=sbr(je)
w=dpa(je)
t=dma(je)
e=exr(je)
60 s1=1.d0-s*s
sp=1.d0+s*u
sm=1.d0-s*u
ij=i+n4*(j-1)
if (ki.ne.1) go to 52
am(ij)=(u1/sp+t)*e
if (je.eq.1) am(ij)=xl
52 if (ki.ne.2) go to 53
am(ij)=u1/sm+w
if (je.eq.1) am(ij)=tau0*xl+u+2.d0*al*uw
53 if (ki.eq.3) am(ij)=(sm+t)*e
if (ki.eq.4) am(ij)=sp+w
if (ki.ne.5) go to 54
am(ij)=u1/sm
if (je.eq.1) am(ij)=1.d0
54 if (ki.ne.6) go to 55
am(ij)=u1*e/sp
if (je.eq.1) am(ij)=-u
55 if (ki.eq.7) am(ij)=sp
if (ki.eq.8) am(ij)=sm*e
if (ki.ne.9) go to 56
am(ij)=t*e
if (je.eq.1) am(ij)=xl
56 if (ki.ne.10) go to 57
am(ij)=w
if (je.eq.1) am(ij)=tau0*xl+u+2.d0*al*uw
57 if (ki.eq.11) am(ij)=(s1/sp+t)*e
if (ki.eq.12) am(ij)=s1/sm+w
if (ki.ne.13) go to 58
am(ij)=0.d0
if (je.eq.1) am(ij)=1.d0
58 if (ki.ne.14) go to 59
am(ij)=0.d0
if (je.eq.1) am(ij)=-u
59 if (ki.eq.15) am(ij)=s1/sm
if (ki.eq.16) am(ij)=s1/sp*e
50 continue
40 continue
x8=(dmm+.5d0*alb*uinc*(fl+fr))*ext
do 70 i=1,n4
if (i.le.nq) bm(i)=(-(bb(1)*uu(i)*uu(i)+bb(2))/
& (1.d0+uu(i)/uinc))*ext+x8
if (i.gt.nq.and.i.le.n2) bm(i)=-(bb(1)*uu(i-nq)*uu(i-nq)+
& bb(2))/(1.d0-uu(i-nq)/uinc)
if (i.gt.n2.and.i.le.n3) bm(i)=-bb(3)/
& (1.d0+uu(i-n2)/uinc)*ext+x8
if (i.gt.n3) bm(i)=-bb(3)/(1.d0-uu(i-n3)/uinc)
70 continue
71 call linsys(am,bm,n4,mis)
return
end
subroutine coefncon(bm,alamb,tau0,uinc,alb,fl,fr,bb,kp,am,
& zal,val,zbe,vbe,sbr,sbl,uu,cc,nq)
implicit real*8(a-h,o-z)
parameter (nqmax=200)
c this code finds the coefficients a(i),b(i),c(i) and d(i)
c stored in array bm(4*nq)
c for 2*2 Rayleigh scattering
c if the scattering is not conservative
c Viik 87 12 22
c
c
logical l1,l2,l3,l4,l5,l6,l8
dimension jmin(16),jmax(16),imin(16),imax(16)
dimension bb(kp),zal(nq),val(nq),zbe(nq),vbe(nq),
& am(1),bm(1),exl(nqmax),exr(nqmax),dpb(nqmax),dmb(nqmax),
& dpa(nqmax),dma(nqmax),sbr(nq),sbl(nq),uu(nq),cc(nq)
stol=1.d0-1.d-12
l6=alamb.lt.stol
n4=4*nq
n3=3*nq
n2=2*nq
c if (tau0.gt.1.d5) go to 130
do 65 i=1,nq
exl(i)=dexp(-tau0*sbl(i))
65 exr(i)=dexp(-tau0*sbr(i))
ext=dexp(-tau0/uinc)
ij=0
do 99 i=1,4
do 99 j=1,4
ij=ij+1
imin(ij)=(i-1)*nq+1
jmin(ij)=(j-1)*nq+1
imax(ij)=i*nq
99 jmax(ij)=j*nq
if (alb.gt.1.d-12) go to 10
do 100 ki=1,16
j1=jmin(ki)
j2=jmax(ki)
i1=imin(ki)
i2=imax(ki)
l1=ki.eq.1.or.ki.eq.3.or.ki.eq.6.or.ki.eq.8
l2=ki.eq.2.or.ki.eq.4.or.ki.eq.5.or.ki.eq.7
l3=ki.eq.9.or.ki.eq.11.or.ki.eq.14.or.ki.eq.16
l4=ki.eq.10.or.ki.eq.12.or.ki.eq.13.or.ki.eq.15
l5=ki.eq.1.or.ki.eq.5.or.ki.eq.9.or.ki.eq.13
do 101 i=i1,i2
ie=i-i1+1
u=uu(ie)
u2=u*u
do 101 j=j1,j2
je=j-j1+1
if (j.gt.n2) go to 102
s=sbl(je)
e=exl(je)
z=zbe(je)
v=vbe(je)
go to 103
102 s=sbr(je)
e=exr(je)
z=zal(je)
v=val(je)
103 sp=1.d0+u*s
sm=1.d0-u*s
uz=u2+z
ij=i+n4*(j-1)
if (l1) am(ij)=uz*e/sp
if (l2) am(ij)=uz/sm
if (l3) am(ij)=v*e/sp
if (l4) am(ij)=v/sm
if (l6) go to 101
if (l5.and.je.eq.1) am(ij)=1.d0
if (ki.eq.2.and.je.eq.1) am(ij)=tau0+u
if (ki.eq.6.and.je.eq.1) am(ij)=-u
if (ki.eq.10.and.je.eq.1) am(ij)=tau0+u
if (ki.eq.14.and.je.eq.1) am(ij)=-u
101 continue
100 continue
do 104 i=1,n4
if (i.le.nq) bm(i)=-(bb(1)*uu(i)*uu(i)+bb(2))*
& ext/(1.d0+uu(i)/uinc)
if (i.gt.nq.and.i.le.n2) bm(i)=-(bb(1)*uu(i-nq)*uu(i-nq)+
& bb(2))/(1.d0-uu(i-nq)/uinc)
if (i.gt.n2.and.i.le.n3) bm(i)=-bb(3)*ext/
& (1.d0+uu(i-n2)/uinc)
if (i.gt.n3) bm(i)=-bb(3)/(1.d0-uu(i-n3)/uinc)
104 continue
201 call linsys(am,bm,n4,kes)
return
c 130 ij=0
c do 131 i=1,2
c do 131 j=1,2
c ij=ij+1
c imin(ij)=(i-1)*nq+1
c imax(ij)=i*nq
c jmin(ij)=(j-1)*nq+1
c 131 jmax(ij)=j*nq
c do 132 ki=1,4
c j1=jmin(ki)
c j2=jmax(ki)
c i1=imin(ki)
c i2=imax(ki)
c do 133 i=i1,i2
c ie=i-i1+1
c u=uu(ie)
c u2=u*u
c do 133 j=j1,j2
c je=j-j1+1
c ij=i+n2*(j-1)
c if (ki.eq.1) am(ij)=(u2+zbe(je))/(1.d0-u*sbl(je))
c if (ki.eq.2) am(ij)=(u2+zal(je))/(1.d0-u*sbr(je))
c if (ki.eq.3) am(ij)=vbe(je)/(1.d0-u*sbl(je))
c if (ki.eq.4) am(ij)=val(je)/(1.d0-u*sbr(je))
c if (l6) go to 133
c if (ki.eq.1.and.je.eq.1) am(ij)=1.d0
c if (ki.eq.3.and.je.eq.1) am(ij)=1.d0
c 133 continue
c 132 continue
c do 134 i=1,n2
c if (i.le.nq) bm(i)=-(bb(1)*uu(i)*uu(i)+bb(2))/
c & (1.d0-uu(i)/uinc)
c if (i.gt.nq) bm(i)=-bb(3)/(1.d0-uu(i-nq)/uinc)
c 134 continue
c call linsys(am,bm,n2,kes)
c 300 format(1x,5e13.5)
c print 300,(bm(klo),klo=1,2*nq)
c return
10 al=alb
uw=1.d0/3.d0
xl=1.d0-al
do 30 j=1,nq
s1=0.d0
s2=0.d0
s3=0.d0
s4=0.d0
z2=sbl(j)
z3=sbr(j)
do 20 i=1,nq
z5=z2*uu(i)
z4=z3*uu(i)
s1=s1+cc(i)/(1.d0+z5)
s2=s2+cc(i)/(1.d0-z5)
s3=s3+cc(i)/(1.d0+z4)
s4=s4+cc(i)/(1.d0-z4)
20 continue
x2=1.d0
if (z2.ge.1.d-12) x2=1.d0/z2
x3=1.d0/z3
dbp1=(1.d0-s1)*x2
dbp3=x2*(uw-x2*(.5d0-dbp1))
dap1=(1.d0-s3)*x3
dap3=x3*(uw-x3*(.5d0-dap1))
dbm1=(s2-1.d0)*x2
dbm3=x2*(x2*(dbm1-.5d0)-uw)
dam1=(s4-1.d0)*x3
dam3=x3*(x3*(dam1-.5d0)-uw)
dpb(j)=(dbp3+dbp1*(zbe(j)+vbe(j)))*al
dmb(j)=(dbm3+dbm1*(zbe(j)+vbe(j)))*al
dpa(j)=(dap3+dap1*(zal(j)+val(j)))*al
dma(j)=(dam3+dam1*(zal(j)+val(j)))*al
30 continue
s5=0.d0
do 11 i=1,nq
11 s5=s5+cc(i)/(1.d0-uu(i)/uinc)
dmu1=(s5-1.d0)*uinc
dmu3=uinc*(uinc*(dmu1-.5d0)-uw)
dmm=(bb(1)*dmu3+(bb(2)+bb(3))*dmu1)*al
do 40 ki=1,16
j1=jmin(ki)
j2=jmax(ki)
i1=imin(ki)
i2=imax(ki)
do 60 i=i1,i2
ie=i-i1+1
u=uu(ie)
u2=u*u
do 60 j=j1,j2
je=j-j1+1
l8=je.eq.1.and..not.l6
if (j.gt.n2) go to 70
s=sbl(je)
e=exl(je)
z=zbe(je)
v=vbe(je)
w=dpb(je)
t=dmb(je)
go to 80
70 s=sbr(je)
e=exr(je)
z=zal(je)
v=val(je)
w=dpa(je)
t=dma(je)
80 sp=1.d0+u*s
sm=1.d0-u*s
uz=u2+z
ij=i+n4*(j-1)
if (ki.ne.1) go to 81
am(ij)=(uz/sp-t)*e
if (l8) am(ij)=xl
81 if (ki.ne.2) go to 82
am(ij)=uz/sm-w
if (l8) am(ij)=u+tau0*xl+2.d0*al*uw
82 if (ki.eq.3) am(ij)=(uz/sp-t)*e
if (ki.eq.4) am(ij)=uz/sm-w
if (ki.ne.5) go to 83
am(ij)=uz/sm
if (l8) am(ij)=1.d0
83 if (ki.ne.6) go to 84
am(ij)=uz/sp*e
if (l8) am(ij)=-u
84 if (ki.eq.7) am(ij)=uz/sm
if (ki.eq.8) am(ij)=uz/sp*e
if (ki.ne.9) go to 85
am(ij)=(v/sp-t)*e
if (l8) am(ij)=xl
85 if (ki.ne.10) go to 86
am(ij)=v/sm-w
if (l8) am(ij)=u+tau0*xl+2.d0*al*uw
86 if (ki.eq.11) am(ij)=(v/sp-t)*e
if (ki.eq.12) am(ij)=v/sm-w
if (ki.ne.13) go to 87
am(ij)=v/sm
if (l8) am(ij)=1.d0
87 if (ki.ne.14) go to 88
am(ij)=v/sp*e
if (l8) am(ij)=-u
88 if (ki.eq.15) am(ij)= v/sm
if (ki.eq.16) am(ij)=v/sp*e
60 continue
40 continue
x8=(dmm+.5d0*alb*uinc*(fl+fr))*ext
do 50 i=1,n4
if (i.le.nq) bm(i)=(-(bb(1)*uu(i)*uu(i)+bb(2))/
& (1.d0+uu(i)/uinc))*ext+x8
if (i.gt.nq.and.i.le.n2) bm(i)=-(bb(1)*uu(i-nq)*uu(i-nq)+
& bb(2))/(1.d0-uu(i-nq)/uinc)
if (i.gt.n2.and.i.le.n3) bm(i)=(-bb(3)/
& (1.d0+uu(i-n2)/uinc))*ext+x8
if (i.gt.n3) bm(i)=-bb(3)/(1.d0-uu(i-n3)/uinc)
50 continue
call linsys(am,bm,n4,kes)
return
end
subroutine intarblr(aiul,aidl,aiur,aidr,alamb,tau0,qs1,
& ni,us,uinc,alb,fl,fr,tau,bb,kp,bm,zal,val,zbe,vbe,
& sbr,sbl,uu,cc,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
c Viik 87 12 22
logical mf,mg, mh
dimension bb(kp),bm(1),zbe(1),zal(1),vbe(1),
& val(1),aiul(ni),aidl(ni),aiur(ni),aidr(ni),
& us(ni),el1(nqmax),el2(nqmax),el3(nqmax),er1(nqmax),
& er2(nqmax),er3(nqmax),fp51(nqmax),fp52(nqmax),
& fp41(nqmax),fp42(nqmax),fm51(nqmax),fm52(nqmax),
& fm41(nqmax),fm42(nqmax),dpb(nqmax),dpa(nqmax),
& dma(nqmax),dmb(nqmax),
& sbr(nq),sbl(nq),uu(nq),cc(nq)
aks=.375d0*alamb*qs1
dze=.25d0*alamb*(1.d0-qs1)
da=aks+dze
stol=1.d0-1.d-14
mf=alamb.ge.stol.and.qs1.ge.stol
mg=alamb.ge.stol.and.qs1.lt.stol
n2=2*nq
n3=3*nq
un2=uinc*uinc
rofr=.5d0* (aks*(un2*fl+fr)+dze*(fl+fr))
5 ta=tau0-tau
do 140 k=1,nq
bel=sbl(k)
ber=sbr(k)
el1(k)=dexp(-bel*tau)
el2(k)=dexp(-bel*ta)
el3(k)=dexp(-bel*tau0)
er1(k)=dexp(-ber*tau)
er2(k)=dexp(-ber*ta)
er3(k)=dexp(-ber*tau0)
140 continue
et=dexp(-tau/uinc)
eta=dexp(-tau0/uinc)
do 100 j=1,ni
v=us(j)
v2=v*v
wf2=aks*v2+dze
ev=dexp(-tau/v)
etv=dexp(-tau0/v)
etav=dexp(-ta/v)
fp=uinc*(et-eta*etav)/(v+uinc)
fm=uinc*(ev-et)/(v-uinc)
ro0=2.d0*(1.d0-un2)+v2*(3.d0*un2-2.d0)
do 110 k=1,nq
ber=sbr(k)
bel=sbl(k)
vka=v*bel
vk=v*ber
belp=1.d0+vka
belm=1.d0-vka
berp=1.d0+vk
berm=1.d0-vk
fp51(k)=(el1(k)-el3(k)*etav)/belp
fp52(k)=(el2(k)-etav)/belm
fp41(k)=(er1(k)-er3(k)*etav)/berp
fp42(k)=(er2(k)-etav)/berm
fm51(k)=(el1(k)-ev)/belm
fm52(k)=(el2(k)-el3(k)*ev)/belp
fm41(k)=(er1(k)-ev)/berm
fm42(k)=(er2(k)-er3(k)*ev)/berp
110 continue
s1=0.d0
s5=0.d0
s9=0.d0
s13=0.d0
do 120 i=1,nq
u=uu(i)
u2=u*u
ui=u/uinc
rp11=0.d0
rp21=0.d0
rp31=0.d0
rp41=0.d0
roi=2.d0*(1.d0-u2)+v2*(3.d0*u2-2.d0)
wf1=aks*roi+dze
wf3=aks*u2+dze
rm11=0.d0
rm21=0.d0
rm31=0.d0
rm41=0.d0
do 130 k=1,nq
ber=sbr(k)
bel=sbl(k)
uka=u*bel
uk=u*ber
belp=1.d0+uka
belm=1.d0-uka
berp=1.d0+uk
berm=1.d0-uk
if (mf) go to 30
go to 10
30 uq=1.d0-u2
qb=1.d0-ber*ber
q11=bm(k)*uq/belp
if (k.eq.1) q11=0.d0
q12=bm(k+nq)*uq/belm
if (k.eq.1) q12=0.d0
q13=bm(k+n2)*berm
q14=bm(k+n3)*berp
q21=bm(k)*uq/belm
if (k.eq.1) q21=0.d0
q22=bm(k+nq)*uq/belp
if (k.eq.1) q22=0.d0
q23=bm(k+n2)*berp
q24=bm(k+n3)*berm
q31=0.d0
q32=0.d0
q33=bm(k+n2)*qb/berp
q34=bm(k+n3)*qb/berm
q41=0.d0
q42=0.d0
q43=bm(k+n2)*qb/berm
q44=bm(k+n3)*qb/berp
go to 20
10 mh=mg.and.k.eq.1
zkb=u2+zbe(k)
zka=u2+zal(k)
q11=bm(k)*zkb/belp
if (mh) q11=0.d0
q12=bm(k+nq)*zkb/belm
if (mh) q12=0.d0
q13=bm(k+n2)*zka/berp
q14=bm(k+n3)*zka/berm
q21=bm(k)*zkb/belm
if (mh) q21=0.d0
q22=bm(k+nq)*zkb/belp
if (mh) q22=0.d0
q23=bm(k+n2)*zka/berm
q24=bm(k+n3)*zka/berp
q31=bm(k)*vbe(k)/belp
if (mh) q31=0.d0
q32=bm(k+nq)*vbe(k)/belm
if (mh) q32=0.d0
q33=bm(k+n2)*val(k)/berp
q34=bm(k+n3)*val(k)/berm
q41=bm(k)*vbe(k)/belm
if (mh) q41=0.d0
q42=bm(k+nq)*vbe(k)/belp
if (mh) q42=0.d0
q43=bm(k+n2)*val(k)/berm
q44=bm(k+n3)*val(k)/berp
20 rp11=rp11+q11*fp51(k)+q12*fp52(k)+q13*fp41(k)+
& q14*fp42(k)
rp31=rp31+q31*fp51(k)+q32*fp52(k)+q33*fp41(k)+
& q34*fp42(k)
rp21=rp21+q21*fp51(k)+q22*fp52(k)+q23*fp41(k)+
& q24*fp42(k)
rp41=rp41+q41*fp51(k)+q42*fp52(k)+q43*fp41(k)+
& q44*fp42(k)
rm11=rm11+q11*fm51(k)+q12*fm52(k)+q13*fm41(k)+
& q14*fm42(k)
rm21=rm21+q21*fm51(k)+q22*fm52(k)+q23*fm41(k)+
& q24*fm42(k)
rm31=rm31+q31*fm51(k)+q32*fm52(k)+q33*fm41(k)+
& q34*fm42(k)
rm41=rm41+q41*fm51(k)+q42*fm52(k)+q43*fm41(k)+
& q44*fm42(k)
130 continue
ax=bb(1)*u2+bb(2)
q15=ax/(1.d0+ui)
q25=ax/(1.d0-ui)
q35=bb(3)/(1.d0+ui)
q45=bb(3)/(1.d0-ui)
rp11=rp11+q15*fp
rp21=rp21+q25*fp
rp31=rp31+q35*fp
rp41=rp41+q45*fp
rm11=rm11+q15*fm
rm21=rm21+q25*fm
rm31=rm31+q35*fm
rm41=rm41+q45*fm
rp1=rp11+rp21
rp3=rp31+rp41
rm1=rm11+rm21
rm3=rm31+rm41
s1=s1+cc(i)*(wf1*rp1+wf2*rp3)
s5=s5+cc(i)*(wf1*rm1+wf2*rm3)
s9=s9+cc(i)*(da*rp3+wf3*rp1)
s13=s13+cc(i)*(da*rm3+wf3*rm1)
120 continue
rofl=.5d0*(aks*(ro0*fl+v2*fr)+dze*(fl+fr))
aiul(j)=s1+rofl*fp
aidl(j)=s5+rofl*fm
aiur(j)=s9+rofr*fp
aidr(j)=s13+rofr*fm
if (alamb.lt.stol) go to 100
x1=1.d0-etav
x2=tau+v-etav*(tau0+v)
y1=1.d0-ev
y2=tau-v+v*ev
slp=bm(1)*x1+bm(nq+1)*x2
slm=bm(1)*y1+bm(nq+1)*y2
aiul(j)=aiul(j)+slp
aidl(j)=aidl(j)+slm
aiur(j)=aiur(j)+slp
aidr(j)=aidr(j)+slm
100 continue
if (alb.lt.1.d-8) return
al=alb
uw=1.d0/3.d0
do 11 j=1,nq
dbp1=0.d0
dbp3=0.d0
dap1=0.d0
dap3=0.d0
dbm1=0.d0
dbm3=0.d0
dam1=0.d0
dam3=0.d0
do 12 i=1,nq
u=uu(i)
u3=u*u*u
z5=sbl(j)
z4=sbr(j)
z55=z5*u
z44=z4*u
x6=cc(i)*u
x7=cc(i)*u3
zm5=1.d0-z55
zp5=1.d0+z55
zm4=1.d0-z44
zp4=1.d0+z44
dbp1=dbp1+x6/zp5
dbp3=dbp3+x7/zp5
dap1=dap1+x6/zp4
dap3=dap3+x7/zp4
dbm1=dbm1+x6/zm5
dbm3=dbm3+x7/zm5
dam1=dam1+x6/zm4
12 dam3=dam3+x7/zm4
dpb(j)=dbp3+dbp1*(zbe(j)+vbe(j))
dmb(j)=dbm3+dbm1*(zbe(j)+vbe(j))
dpa(j)=dap3+dap1*(zal(j)+val(j))
dma(j)=dam3+dam1*(zal(j)+val(j))
if (.not.mf) go to 11
dmb(j)=dbm1-dbm3
dpb(j)=dbp1-dbp3
dma(j)=.5d0+uw*z4+(1.d0-z4*z4)*dam1
dpa(j)=.5d0-uw*z4+(1.d0-z4*z4)*dap1
11 continue
s5=0.d0
do 13 i=1,nq
13 s5=s5+cc(i)/(1.d0-uu(i)/uinc)
dmu1=(s5-1.d0)*uinc
dmu3=uinc*(uinc*(dmu1-.5d0)-uw)
dmm=bb(1)*dmu3+(bb(2)+bb(3))*dmu1
mo=1
if (alamb.gt.1.d0-1.d-10) mo=2
s1=0.d0
s2=0.d0
do 14 i=mo,nq
s1=s1+bm(i)*dmb(i)*el3(i)
14 s2=s2+bm(i+nq)*dpb(i)
s3=0.d0
s4=0.d0
do 15 i=1,nq
s3=s3+bm(i+n2)*dma(i)*er3(i)
15 s4=s4+bm(i+n3)*dpa(i)
ag=0.d0
if (alamb.gt.1.d0-1.d-10)
& ag=bm(1)+bm(nq+1)*(tau0-2.d0*uw)
af=al*(ag+s1+s2+s3+s4+dmm*eta)+
& .5d0*alb*uinc*eta*(fl+fr)
do 16 i=1,ni
ex=dexp(-ta/us(i))
aiul(i)=aiul(i)+af*ex
16 aiur(i)=aiur(i)+af*ex
return
end
subroutine intgausscon(aiul,aidl,aiur,aidr,tau0,uinc,
& tau,bb,kp,bm,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 scattering in homogeneous
c plane-parallel conservative atmosphere
c Viik 87 12 23
dimension bb(kp),bm(1),aiul(1),aidl(1),
& aiur(1),aidr(1),dpa(nqmax),dma(nqmax),dpb(nqmax),dmb(nqmax),
& e51(nqmax),e52(nqmax),e41(nqmax),e42(nqmax),
& sbr(nq),sbl(nq),uu(nq)
n2=2*nq
n3=3*nq
5 ta=tau0-tau
ex=dexp(-tau/uinc)
do 130 k=1,nq
bl=sbl(k)
br=sbr(k)
e51(k)=dexp(-bl*tau)
e52(k)=dexp(-bl*ta)
e41(k)=dexp(-br*tau)
130 e42(k)=dexp(-br*ta)
do 100 i=1,nq
v=uu(i)
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
if (tau0.gt.1.d5) go to 170
do 110 j=2,nq
bl=sbl(j)
e1=e51(j)
e2=e52(j)
s1=s1+bm(j)*e1/(1.d0+v*bl)
s2=s2+bm(nq+j)*e2/(1.d0-v*bl)
s5=s5+bm(j)*e1/(1.d0-v*bl)
110 s6=s6+bm(nq+j)*e2/(1.d0+v*bl)
do 120 j=1,nq
br=sbr(j)
e1=e41(j)
e2=e42(j)
dp=1.d0+br*v
dm=1.d0-br*v
dx=1.d0-br*br
s3=s3+bm(n2+j)*dm*e1
s4=s4+bm(n3+j)*dp*e2
s7=s7+bm(n2+j)*dp*e1
s8=s8+bm(n3+j)*dm*e2
s9=s9+bm(n2+j)*e1/dp*dx
s10=s10+bm(n3+j)*e2/dm*dx
s11=s11+bm(n2+j)*e1/dm*dx
120 s12=s12+bm(n3+j)*e2/dp*dx
bq1=bm(nq+1)
go to 200
170 do 180 j=2,nq
bl=sbl(j)
e1=e51(j)
s1=s1+bm(j)*e1/(1.d0+v*bl)
180 s5=s5+bm(j)*e1/(1.d0-v*bl)
do 190 j=1,nq
br=sbr(j)
e1=e41(j)
dp=1.d0+br*v
dm=1.d0-br*v
dx=1.d0-br*br
s3=s3+bm(j+nq)*dm*e1
s7=s7+bm(j+nq)*dp*e1
s9=s9+bm(j+nq)*e1/dp*dx
190 s11=s11+bm(j+nq)*e1/dm*dx
bq1=0.d0
200 v1=1.d0-v*v
rx=bb(1)*v*v+bb(2)
ax=rx/(1.d0+v/uinc)
bx=rx/(1.d0-v/uinc)
cx=bb(3)/(1.d0+v/uinc)
dx=bb(3)/(1.d0-v/uinc)
aiul(i)=bm(1)+v1*s1+bq1*(tau+v)+v1*s2+s3+s4+ax*ex
aidl(i)=bm(1)+v1*s5+bq1*(tau-v)+v1*s6+s7+s8+bx*ex
aiur(i)=bm(1)+bq1*(tau+v)+s9+s10+cx*ex
aidr(i)=bm(1)+bq1*(tau-v)+s11+s12+dx*ex
100 continue
return
end
subroutine intgaussncon(aiul,aidl,aiur,aidr,alamb,tau0,
& uinc,tau,bb,kp,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 scattering
c if scattering is not conservative
c Viik 87 12 22
logical l1
dimension bb(kp),bm(1),zbe(1),zal(1),vbe(1),
& val(1),aiul(1),aidl(1),aiur(1),aidr(1),
& e51(nqmax),e52(nqmax),e41(nqmax),e42(nqmax),
& sbr(nq),sbl(nq),uu(nq)
stol=1.d0-1.d-12
l1=alamb.ge.stol
rg1=0.d0
rg2=0.d0
if (.not.l1) go to 50
rg1=bm(1)
rg2=bm(nq+1)
50 n2=2*nq
n3=3*nq
5 ta=tau0-tau
ex=dexp(-tau/uinc)
do 130 k=1,nq
bl=sbl(k)
br=sbr(k)
e51(k)=dexp(-bl*tau)
e52(k)=dexp(-bl*ta)
e41(k)=dexp(-br*tau)
130 e42(k)=dexp(-br*ta)
do 100 i=1,nq
v=uu(i)
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
s13=0.d0
s14=0.d0
s15=0.d0
s16=0.d0
if (tau0.gt.1.d5) go to 170
do 110 j=1,nq
bl=sbl(j)
el1=e51(j)
el2=e52(j)
br=sbr(j)
er1=e41(j)
er2=e42(j)
ql=v*v+zbe(j)
qr=v*v+zal(j)
vbl=v*bl
vbr=v*br
bml=1.d0-vbl
bpl=1.d0+vbl
bmr=1.d0-vbr
bpr=1.d0+vbr
if (l1.and.j.eq.1) go to 20
fs1=bm(j)*ql*el1/bpl
fs2=bm(j+nq)*ql*el2/bml
fs5=bm(j)*ql*el1/bml
fs6=bm(j+nq)*ql*el2/bpl
fs9=bm(j)*vbe(j)*el1/bpl
fs10=bm(j+nq)*vbe(j)*el2/bml
fs13=bm(j)*vbe(j)*el1/bml
fs14=bm(j+nq)*vbe(j)*el2/bpl
go to 10
20 fs1=0.d0
fs2=0.d0
fs5=0.d0
fs6=0.d0
fs9=0.d0
fs10=0.d0
fs13=0.d0
fs14=0.d0
10 s1=s1+fs1
s2=s2+fs2
s3=s3+bm(j+n2)*qr*er1/bpr
s4=s4+bm(j+n3)*qr*er2/bmr
s5=s5+fs5
s6=s6+fs6
s7=s7+bm(j+n2)*qr*er1/bmr
s8=s8+bm(j+n3)*qr*er2/bpr
s9=s9+fs9
s10=s10+fs10
s11=s11+bm(j+n2)*val(j)*er1/bpr
s12=s12+bm(j+n3)*val(j)*er2/bmr
s13=s13+fs13
s14=s14+fs14
s15=s15+bm(j+n2)*val(j)*er1/bmr
s16=s16+bm(j+n3)*val(j)*er2/bpr
110 continue
go to 200
170 do 180 j=1,nq
bl=sbl(j)
br=sbr(j)
el=e51(j)
er=e41(j)
ql=v*v+zbe(j)
qr=v*v+zal(j)
vbl=v*bl
vbr=v*br
bml=1.d0-vbl
bpl=1.d0+vbl
bmr=1.d0-vbr
bpr=1.d0+vbr
if (l1.and.j.eq.1) go to 30
fs1=bm(j)*ql*el/bpl
fs5=bm(j)*ql*el/bml
fs9=bm(j)*vbe(j)*el/bpl
fs13=bm(j)*vbe(j)*el/bml
go to 40
30 fs1=0.d0
fs5=0.d0
fs9=0.d0
fs13=0.d0
40 s1=s1+fs1
s3=s3+bm(j+nq)*qr*er/bpr
s5=s5+fs5
s7=s7+bm(j+nq)*qr*er/bmr
s9=s9+fs9
s11=s11+bm(j+nq)*val(j)*er/bpr
s13=s13+fs13
180 s15=s15+bm(j+nq)*val(j)*er/bmr
rg2=0.d0
200 rx=bb(1)*v*v+bb(2)
ax=rx/(1.d0+v/uinc)
bx=rx/(1.d0-v/uinc)
cx=bb(3)/(1.d0+v/uinc)
dx=bb(3)/(1.d0-v/uinc)
aiul(i)=rg1+s1+rg2*(tau+v)+s2+s3+s4+ax*ex
aidl(i)=rg1+s5+rg2*(tau-v)+s6+s7+s8+bx*ex
aiur(i)=rg1+s9+rg2*(tau+v)+s10+s11+s12+cx*ex
aidr(i)=rg1+s13+rg2*(tau-v)+s14+s15+s16+dx*ex
100 continue
return
end
subroutine gaussint(aium,aqum,auum,avum,pu,
& aidm,aqdm,audm,avdm,pd,aul0,aur0,adl0,adr0,
& am,bm,alamb,qs1,tau0,tau,
& uinc,alb,fi,kfi,fl,fr,fu,fv,zal,val,zbe,vbe,s,a,b,uu,cc,nq)
implicit real*8(a-h,o-z)
parameter(nqmax=200)
c this code finds the vectors of intensity at
c optical depth tau for the Rayleigh
c scattering phase matrix in a homogeneous plane-parallel
c atmosphere
c kfi=1 for the first call
c kfi=2 next calls, if lambda and tau0
c are the same
c kfi=4 next calls, different fi only
c kfi=5 next calls, if lambda is the same
c but tau0 is different
c Viik 87 12 22
logical m1
dimension au(nqmax),ad(nqmax),auu(4,nqmax),add(4,nqmax),
& aiul(nqmax),aidl(nqmax),aiur(nqmax),aidr(nqmax),uf(nqmax),
& aium(nq),aqum(nq),auum(nq),avum(nq),pu(nq),
& aidm(nq),aqdm(nq),audm(nq),avdm(nq),pd(nq),
& aa(9),bb(3),
& am(1),bm(1),zal(nq),val(nq),zbe(nq),vbe(nq),
& sbr(nq),sbl(nq),uu(nq),cc(nq),a(nq),b(nq),s(nq),
& az(4,nqmax),bz(4,nqmax),betaz(4,nqmax),ala(4),
& aul0(nq),aur0(nq),adl0(nq),adr0(nq)
stol=1.d0-1.d-14
tol=1d-12
nmax=200
g=(1.d0-qs1)/(1.d0+2.d0*qs1)
m1=alamb.ge.stol.and.qs1.ge.stol
if (kfi.eq.2) go to 20
if (kfi.eq.4) go to 50
if (kfi.eq.5) go to 70
call zeropol(sbr,sbl,alamb,qs1,uu,cc,nq)
nl=nq
ng=nq
do 93 ia=1,nq
93 uf(ia)=uu(ia)
pi=3.141592653589793d0
do 10 liw=1,4
miw=liw+5
if (liw.eq.1.or.liw.eq.2) ala(liw)=alamb*qs1
if (liw.eq.3.or.liw.eq.4) ala(liw)=alamb*(1.d0-3.d0*g)/
& (1.d0+2.d0*g)
alam=ala(liw)
call zeronat(s,alam,uu,cc,nq,miw)
do 10 i=1,nq
betaz(liw,i)=s(i)
10 continue
70 do 90 liw=1,4
miw=liw+5
do 91 i=1,nq
91 s(i)=betaz(liw,i)
if (tau0.lt.1.d5) call ficoeff(a,b,s,uu,alamb,tau0,nq)
if (tau0.gt.1.d5) call sicoeff(a,s,uu,nq)
do 92 i=1,nq
az(liw,i)=a(i)
92 bz(liw,i)=b(i)
90 continue
20 do 30 liw=1,4
miw=liw+5
alam=ala(liw)
do 40 i=1,nq
s(i)=betaz(liw,i)
a(i)=az(liw,i)
40 b(i)=bz(liw,i)
if (tau0.le.1.d5)
& call fintst(au,ad,alam,tau,tau0,uinc,uf,s,a,b,nq,nq)
if (tau0.gt.1.d5)
& call sintst(au,ad,alam,tau,uinc,uf,s,a,nq,nq)
do 11 i=1,nq
auu(liw,i)=au(i)*4.d0
11 add(liw,i)=ad(i)*4.d0
30 continue
c
c liw=1 corresponds to fii super 1
c liw=2 corresponds to fii super 2
c liw=3 corresponds to v super 0
c liw=4 corresponds to v super 1
c
kp=3
kp2=kp*kp
call albega(zal,val,zbe,vbe,alamb,uinc,qs1,fl,fr,
& aa,bb,kp2,kp,sbr,sbl,uu,cc,nq)
if (m1) call coefcon(bm,tau0,uinc,alb,fl,fr,bb,kp,am,
& sbr,sbl,uu,cc,nq)
if (.not.m1) call coefncon(bm,alamb,tau0,uinc,alb,fl,fr,bb,kp,am,
& zal,val,zbe,vbe,sbr,sbl,uu,cc,nq)
80 if (m1) call intgausscon(aiul,aidl,aiur,aidr,tau0,uinc,
& tau,bb,kp,bm,sbr,sbl,uu,nq)
if (.not.m1) call intgaussncon(aiul,aidl,aiur,aidr,alamb,tau0,
& uinc,tau,bb,kp,bm,zal,val,zbe,
& vbe,sbr,sbl,uu,nq)
do 13 i=1,nq
aul0(i)=aiul(i)
aur0(i)=aiur(i)
adl0(i)=aidl(i)
13 adr0(i)=aidr(i)
50 ui2=uinc*uinc
sinc=dsqrt(1.d0-ui2)
s1=dsin(fi)
s2=dsin(2.d0*fi)
c1=dcos(fi)
c2=dcos(2.d0*fi)
do 60 i=1,nq
v=uu(i)
v2=v*v
sv=dsqrt(1.d0-v2)
vu=v*uinc
c
aul=aiul(i)+sinc*sv*.375d0*auu(1,i)*
& (-2.d0*vu*c1*fl+v*s1*fu)+
& .1875d0*auu(2,i)*(vu*vu*c2*fl-
& v2*c2*fr-v*vu*s2*fu)
c
adl=aidl(i)+sinc*sv*.375d0*add(1,i)*
& (2.d0*vu*c1*fl-v*s1*fu)+
& .1875d0*add(2,i)*(vu*vu*c2*fl-
& v2*c2*fr-v*vu*s2*fu)
c
aur=aiur(i)+.1875d0*auu(2,i)*
& (-ui2*c2*fl+c2*fr+uinc*s2*fu)
c
adr=aidr(i)+.1875d0*add(2,i)*
& (-ui2*c2*fl+c2*fr+uinc*s2*fu)
c
aium(i)=aul+aur
aidm(i)=adl+adr
aqum(i)=aur-aul
aqdm(i)=adr-adl
c
auum(i)=.375d0*sinc*sv*auu(1,i)*
& (2.d0*uinc*s1*fl+c1*fu)+
& .1875d0*auu(2,i)*(-2.d0*v*ui2*s2*fl+
& 2.d0*v*s2*fr-2.d0*vu*c2*fu)
c
audm(i)=.375d0*sinc*sv*add(1,i)*
& (2.d0*uinc*s1*fl+c1*fu)+
& .1875d0*add(2,i)*(2.d0*v*ui2*s2*fl-
& 2.d0*v*s2*fr+2.d0*vu*c2*fu)
c
avum(i)=.375d0*fv*(-vu*auu(3,i)+
& sinc*sv*auu(4,i)*c1)
c
avdm(i)=.375d0*fv*(vu*add(3,i)+
& sinc*sv*add(4,i)*c1)
c
anim=aium(i)
alug=aqum(i)
pu(i)=10.d0
if (anim.ge.tol) pu(i)=dsqrt(aqum(i)*aqum(i)+
& auum(i)*auum(i)+avum(i)*avum(i))/anim
anim=aidm(i)
alug=aqdm(i)
pd(i)=10.d0
if (anim.ge.tol) pd(i)=dsqrt(aqdm(i)*aqdm(i)+
& audm(i)*audm(i)+avdm(i)*avdm(i))/anim
60 continue
return
end
subroutine stokes(aium,aqum,auum,avum,pu,
& aidm,aqdm,audm,avdm,pd,
& aul0,aur0,adl0,adr0,
& tau0,ni,us,kinc,inc,fl,fr,fu,fv,
& fi,alb,uu,cc,s,a,b,nq)
implicit real*8(a-h,o-z)
parameter(nqmax=200)
c
c subroutine STOKES finds the vectors of
c emergent intensities at tau=0 and tau=tau0
c for the RAYLEIGH scattering in a homogeneous
c plane-parallel atmosphere (lambda=1).
c
dimension us(ni),aium(ni),aqum(ni),auum(ni),avum(ni),pu(ni),
& aidm(ni),aqdm(ni),audm(ni),avdm(ni),pd(ni),
& xt(5,nqmax),yt(5,nqmax),psiit(nqmax),aksiit(nqmax),
& fiit(nqmax),etat(nqmax),hiit(nqmax),sigt(nqmax),dzet(nqmax),
& tett(nqmax),psii(nqmax),aksii(nqmax),fii(nqmax),eta(nqmax),
& hii(nqmax),sig(nqmax),dze(nqmax),tet(nqmax),
& xf(nqmax),yf(nqmax),xg(5,nqmax),yg(5,nqmax),
& xp(5,nqmax),yp(5,nqmax),w1(nqmax),w2(nqmax),wm1(nqmax),
& wm2(nqmax),
& wm3(nqmax),wm4(nqmax),gl1(nqmax),gr1(nqmax),vsl(nqmax),
& vsr(nqmax),
& w6(nqmax),w7(nqmax),w8(nqmax),w9(nqmax),
& uu(nq),cc(nq),s(nq),a(nq),b(nq),
& aul0(ni),aur0(ni),adl0(ni),adr0(ni)
alamb=1.d0
tol=1d-12
ing=1
if (kinc.ne.1) go to 260
pi=3.14159265358979d0
do 10 liw=1,5
miw=liw+5
call zeronat(s,alamb,uu,cc,nq,miw)
call ficoeff(a,b,s,uu,alamb,tau0,nq)
do 16 i=1,nq
u=uu(i)
call xyfun(x,y,u,tau0,s,a,b,nq,ing)
xf(i)=x
yf(i)=y
16 continue
do 15 i=1,nq
xg(liw,i)=xf(i)
15 yg(liw,i)=yf(i)
do 25 i=1,ni
u18=us(i)
call xyfun(x8,y8,u18,tau0,s,a,b,nq,ing)
call dxidyi(dxu,dyu1,tau0,0d0,u18,s,a,b,nq)
call dxidyi(dxu1,dyu,tau0,tau0,u18,s,a,b,nq)
xt(liw,i)=dxu
yt(liw,i)=dyu
xp(liw,i)=x8
25 yp(liw,i)=y8
10 continue
a0=0.d0
a1=0.d0
a2=0.d0
a3=0.d0
b0=0.d0
b1=0.d0
b2=0.d0
b3=0.d0
do 35 i=1,nq
c=cc(i)
u=uu(i)
a0=a0+c*xg(5,i)
a1=a1+c*u*xg(5,i)
a2=a2+c*u*u*xg(5,i)
a3=a3+c*u*u*u*xg(5,i)
b0=b0+c*yg(5,i)
b1=b1+c*u*yg(5,i)
b3=b3+c*u*u*u*yg(5,i)
35 b2=b2+c*u*u*yg(5,i)
cu=4.d0/3.d0/(a0-a2)
sq7=(b0-b2)/(a1-a3+b1-b3)
sq8=(4.d0/3.d0-a0+a2)/(a1+b1-a3-b3)
do 45 i=1,ni
u=us(i)
cfg=sq7*(xp(5,i)+yp(5,i))
cfgt=sq7*(xt(5,i)+yt(5,i))
xt(5,i)=xt(5,i)+cfg+u*cfgt
xp(5,i)=xp(5,i)+cfg*u
yt(5,i)=yt(5,i)-cfg-cfgt*u
45 yp(5,i)=yp(5,i)-cfg*u
do 240 i=1,nq
cfg=sq7*(xg(5,i)+yg(5,i))*uu(i)
xg(5,i)=xg(5,i)+cfg
240 yg(5,i)=yg(5,i)-cfg
a0=0.d0
a2=0.d0
b0=0.d0
b2=0.d0
do 250 i=1,nq
c=cc(i)
u=uu(i)
a0=a0+c*xg(5,i)
a2=a2+c*u*u*xg(5,i)
b0=b0+c*yg(5,i)
250 b2=b2+c*u*u*yg(5,i)
s1=a0-a2-4.d0/3.d0
s2=b0-b2
55 sa0=0.d0
sb0=0.d0
sa1=0.d0
sb1=0.d0
sa2=0.d0
sb2=0.d0
sa3=0.d0
sb3=0.d0
a1=0.d0
b1=0.d0
a2=0.d0
b2=0.d0
do 20 i=1,nq
u=uu(i)
c=cc(i)
cu=c*u
cuu=cu*u
cuuu=cuu*u
sa0=sa0+c*xg(4,i)
sb0=sb0+c*yg(4,i)
sa1=sa1+cu*xg(4,i)
sb1=sb1+cu*yg(4,i)
sa2=sa2+cuu*xg(4,i)
sb2=sb2+cuu*yg(4,i)
sa3=sa3+cuuu*xg(4,i)
sb3=sb3+cuuu*yg(4,i)
a1=a1+cu*xg(5,i)
b1=b1+cu*yg(5,i)
a2=a2+cuu*xg(5,i)
20 b2=b2+cuu*yg(5,i)
c0=sa0+sb0-8.d0/3.d0
d0=sa0-sb0-8.d0/3.d0
c1=sa1+sb1
c2=sa2+sb2
d1=sa1-sb1
d2=sa2-sb2
d3=sa3-sb3
ak1=a1+b1
ak2=a2+b2
de1=a1-b1
de2=a2-b2
qk=(c0-c2)/((d0-d2)*tau0+2.d0*(d1-d3))
sd1=1.d0/(d0*de1-d1*de2)
sd2=1.d0/(c0*ak1-c1*ak2-2.d0*qk*
& (d1*ak1-d0*ak2))
g1=2.d0*sd1*(ak1*de1-ak2*de2)
g2=2.d0*sd2*(ak1*de1-ak2*de2)
g3=sd1*(d1*ak1-d0*ak2)
g4=sd2*(c1*de1-c0*de2-2.d0*qk*
& (d0*de1-d1*de2))
g5=sd1*(c1*de1-c0*de2)
g6=sd2*(d1*ak1-d0*ak2)
an1=(g1-g2)/2.d0
an2=(g1+g2)/2.d0
an3=(g3-g4)/2.d0
an4=(g3+g4)/2.d0
u3=(g5-g6)/2.d0
u4=(g5+g6)/2.d0
u5=sd2*(c0*ak1-c1*ak2)
qku=qk*(u4-u3)
uu5=1.d0+2.d0*qku
qka=qk*(an2-an1)
qn2=3.d0/8.d0*qk*(d0-d2)
qn1=qn2*(an2-an1)
ska=1.d0-qn2*((an2-an1)*ak1+
& (u4-u3)*c1-u5*d2)
do 30 i=1,ni
u=us(i)
u2=u*u
xy=xp(4,i)-yp(4,i)
xyt=xt(4,i)-yt(4,i)
psii(i)=u*(an1*yp(5,i)-an2*xp(5,i))
aksii(i)=u*(an2*yp(5,i)-an1*xp(5,i))
fii(i)=(1.d0+an4*u)*xp(5,i)-an3*u*yp(5,i)
eta(i)=(1.d0-an4*u)*yp(5,i)+an3*u*xp(5,i)
hii(i)=(1.d0-u4*u)*xp(4,i)+u3*u*yp(4,i)+
& qku*u2*xy
sig(i)=(1.d0+u4*u)*yp(4,i)-u3*u*xp(4,i)-
& qku*u2*xy
dze(i)=.5d0*u*(an1*yp(4,i)-an2*xp(4,i))+
& .5d0*qka*u2*xy
tet(i)=.5d0*u*(an2*yp(4,i)-an1*xp(4,i))-
& .5d0*qka*u2*xy
psiit(i)=psii(i)/u+u*(an1*yt(5,i)-an2*xt(5,i))
aksiit(i)=aksii(i)/u+u*(an2*yt(5,i)-an1*xt(5,i))
fiit(i)=an4*xp(5,i)+(1.d0+an4*u)*xt(5,i)-
& an3*yp(5,i)-an3*u*yt(5,i)
etat(i)=-an4*yp(5,i)+(1.d0-an4*u)*yt(5,i)+
& an3*xp(5,i)+an3*u*xt(5,i)
hiit(i)=-u4*xp(4,i)+(1.d0-u4*u)*xt(4,i)+
& u3*yp(4,i)+u3*u*yt(4,i)+qku*2.d0*u*xy+
& qku*u2*xyt
sigt(i)=u4*yp(4,i)+(1.d0+u4*u)*yt(4,i)-
& u3*xp(4,i)-u3*u*xt(4,i)-qku*2.d0*u*xy-
& qku*u2*xyt
dzet(i)=.5d0*(an1*yp(4,i)-an2*xp(4,i))+
& .5d0*u*(an1*yt(4,i)-an2*xt(4,i))+
& qka*u*xy+.5d0*qka*u2*xyt
tett(i)=.5d0*(an2*yp( 4,i)-an1*xp(4,i))+
& .5d0*u*(an2*yt(4,i)-an1*xt(4,i))-
& qka*u*xy-.5d0*qka*u2*xyt
gl1(i)=qn1*(xp(5,i)+yp(5,i))
gr1(i)=qn2*((u4-u3)*(xp(4,i)+yp(4,i))-
& u5*u*(xp(4,i)-yp(4,i)))
vsl(i)=1.d0-gl1(i)
vsr(i)=1.d0-gr1(i)
30 continue
260 uinc=us(inc)
ab1=0.5d0*alb*uinc/(1.d0-ska*alb)
do 230 i=1,ni
u=us(i)
w6(i)=ab1*(gl1(i)*gl1(inc)*fl+gl1(i)*gr1(inc)*fr)
w7(i)=ab1*(gr1(i)*gl1(inc)*fl+gr1(i)*gr1(inc)*fr)
w8(i)=ab1*(vsl(i)*gl1(inc)*fl+vsl(i)*gr1(inc)*fr)
230 w9(i)=ab1*(vsr(i)*gl1(inc)*fl+vsr(i)*gr1(inc)*fr)
c=3.d0/32.d0*uinc
sinc=dsqrt(1.d0-uinc*uinc)
s1=dsin(fi)
s2=dsin(2.d0*fi)
c1=dcos(fi)
c2=dcos(2.d0*fi)
do 60 i=1,ni
u=us(i)
su=dsqrt(1.d0-u*u)
uin=uinc+u
uim=u-uinc
aul0(i)=2.d0*c/uin*((psii(i)*psii(inc)+
& 2.d0*fii(i)*fii(inc))*fl+
& (psii(i)*hii(inc)+2.d0*fii(i)*dze(inc))*fr-
& (aksii(i)*aksii(inc)+2.d0*eta(i)*eta(inc))*fl-
& (aksii(i)*sig(inc)+2.d0*eta(i)*tet(inc))*fr)
aur0(i)=2.d0*c/uin*((hii(i)*psii(inc)+
& 2.d0*dze(i)*fii(inc))*fl+
& (hii(i)*hii(inc)+2.d0*dze(i)*dze(inc))*fr-
& (sig(i)*aksii(inc)+2 .d0*tet(i)*eta(inc))*fl-
& (sig(i)*sig(inc)+2.d0*tet(i)*tet(inc))*fr)
w1(i)=xp(1,inc)*yp(1,i)-yp(1,inc)*xp(1,i)
w2(i)=xp(2,inc)*yp(2,i)-yp(2,inc)*xp(2,i)
wm1(i)=xp(1,inc)*xp(1,i)-yp(1,inc)*yp(1,i)
wm2(i)=xp(2,inc)*xp(2,i)-yp(2,inc)*yp(2,i)
wm3(i)=xp(3,i)*xp(3,inc)-yp(3,i)*yp(3,inc)
wm4(i)=xp(4,i)*xp(4,inc)-yp(4,inc)*yp(4,i)
if (u.eq.uinc) go to 200
adl0(i)=2.d0*c/uim*((aksii(i)*psii(inc)+
& 2.d0*eta(i)*fii(inc))*
& fl+(aksii(i)*hii(inc)+2.d0*eta(i)*dze(inc))*fr-
& (psii(i)*aksii(inc)+2.d0*fii(i)*eta(inc))*fl-
& (psii(i)*sig(inc)+2.d0*fii(i)*tet(inc))*fr)
adr0(i)=2.d0*c/uim*((sig(i)*psii(inc)+
& 2.d0*tet(i)*fii(inc))*fl+
& (sig(i)*hii(inc)+2.d0*tet(i)*dze(inc))*fr-
& (hii(i)*aksii(inc)+2.d0*dze(i)*eta(inc))*fl-
& (hii(i)*sig(inc)+2.d0*dze(i)*tet(inc))*fr)
go to 60
200 adl0(i)=2.d0*c*((aksiit(i)*psii(inc)+
& 2.d0*etat(i)*fii(inc))*fl+
& (aksiit(i)*hii(inc)+2.d0*etat(i)*dze(inc))*fr-
& (psiit(i)*aksii(inc)+2.d0*fiit(i)*eta(inc))*fl-
& (psiit(i)*sig(inc)+2.d0*fiit(i)*tet(inc))*fr)
adr0(i)=2.d0*c*((sigt(i)*psii(inc)+2.d0*tett(i)*fii(inc))*fl+
& (sigt(i)*hii(inc)+2.d0*tett(i)*dze(inc))*fr-
& (hiit(i)*aksii( inc)+2.d0*dzet(i)*eta(inc))*fl-
& (hiit(i)*sig(inc)+2.d0*dzet(i)*tet(inc))*fr)
60 continue
do 50 i=1,ni
u=us(i)
uin=uinc+u
uim=u-uinc
uc=u*uinc
su=dsqrt(1.d0-u*u)
r1=4.d0*c*su*sinc*(-2.d0*uc*c1*fl+u*s1*fu)
r2=2.d0*c*(uc*uc*c2*fl-u*u*c2*fr-u*uc*s2*fu)
aul0(i)=aul0(i)+w6(i)
ail=aul0(i)+(r1*wm1(i)+r2*wm2(i))/uin
r3=2.d0*c*(-uinc*uinc*c2*fl+c2*fr+uinc*s2*fu)
aur0(i)=aur0(i)+w7(i)
air=aur0(i)+r3*wm2(i)/uin
r4=4.d0*c*su*sinc*(2.d0*uinc*s1*fl+c1*fu)
r5=4.d0*c*(-uc*uinc*s2*fl+u*s2*fr-uc*c2*fu)
auu=(r4*wm1(i)+r5*wm2(i))/uin
auum(i)=auu
aiu=ail+air
aium(i)=aiu
aqu=air-ail
aqum(i)=aqu
c
avu=4.d0*c/uin*(-uc*wm3(i)+su*sinc*
& wm4(i)*c1)*fv
avum(i)=avu
if (aium(i).lt.tol) then
pu(i)=1d1
else
pu(i)=dsqrt(aqum(i)*aqum(i)+auum(i)*auum(i)+
& avum(i)*avum(i))/aium(i)
endif
if (u.eq.uinc) go to 210
adl0(i)=adl0(i)+w8(i)
aild=adl0(i)-(r1*w1(i)-r2*w2(i))/uim
adr0(i)=adr0(i)+w9(i)
aird=adr0(i)+r3*w2(i)/uim
aud=(r4*w1(i)-r5*w2(i))/uim
audm(i)=aud
aid=aild+aird
aidm(i)=aid
aqd=aird-aild
aqdm(i)=aqd
c
avd=4.d0*c/uim*(uc*(yp(3,i)*xp(3,inc)-xp(3,i)*yp(3,inc))+
& su*sinc*(yp(4,i)*xp(4,inc)-xp(4,i)*yp(4,inc))*c1)*fv
avdm(i)=avd
go to 220
210 adl0(i)=adl0(i)+w8(i)
rt1=xp(1,inc)*yt(1,i)-yp(1,inc)*xt(1,i)
rt2=xp(2,inc)*yt(2,i)-yp(2,inc)*xt(2,i)
aild=adl0(i)-r1*rt1+r2*rt2
adr0(i)=adr0(i)+w9(i)
aird=adr0(i)+r3*rt2
aud=r4*rt1-r5*rt2
audm(i)=aud
aid=aild+aird
aidm(i)=aid
aqd=aird-aild
aqdm(i)=aqd
c
avd=4.d0*c*fv*(uc*(yt(3,i)*xp(3,inc)-xt(3,i)*yp(3,inc))+
& su*sinc*(yt(4,i)*xp(4,inc)-xt(4,i)*yp(4,inc)))
avdm(i)=avd
220 if (aidm(i).lt.tol) then
pd(i)=1d1
else
pd(i)=dsqrt(aqdm(i)*aqdm(i)+audm(i)*audm(i)+
& avdm(i)*avdm(i))/aidm(i)
endif
50 continue
75 continue
5 continue
return
end
subroutine allint(aium,aqum,auum,avum,pu,
& aidm,aqdm,audm,avdm,pd,
& aul0,aur0,adl0,adr0,
& sbr,sbl,am,bm,alamb,tau0,tau,qs1,
& uinc,alb,fi,kfi,fl,fr,fu,fv,a,b,s,uu,cc,us,ni,nq)
implicit real*8(a-h,o-z)
parameter (nqmax=200)
c this code finds all four components of the STOKES vector of
c intensity in (I,Q,U,V) representation at optical depth tau for
c planetary problem in an homogeneous plane-parallel atmosphere.
c Rayleigh-Cabannes scattering.
c Intensities of the upward radiation: aium,aqum,auum,avum;
c downward radiation: aidm,aqdm,audm,avdm;
c pu and pd - degrees of polarization, for radiation moving
c up and down, respectively;
c alamb - albedo of single scattering;
c tau0 - optical thickness of the atmosphere;
c qs1 - depolarization index (qs1=0 - isotropic scattering,
c qs1=1 - Rayleigh scattering);
c uinc - arc cos uinc is the angle of incident beam;
c alb - albedo of the lower boundary of the atmosphere
c (Lambert reflector is assumed);
c fi - angle of azimuth;
c fl, fr, fu, fv - components of the incident flux
c (NB! (l,r) representation);
c uu - array of the Gaussian points;
c cc - array of the Gaussian weigths;
c nq - order of the Gaussian quadrature;
c us - array of arbitrary angular points
c where we want to find the Stokes vector;
c ni - dimension of the us array;
c sbr,sbl,am,bm,a,b,s - auxialiary arrays.
c
c Viik 87 12 22
300 format(4x,5d15.6)
logical m1
dimension au(nqmax),ad(nqmax),auu(4*nqmax),add(4*nqmax),
& aiul(nqmax),aidl(nqmax),aiur(nqmax),aidr(nqmax),us(ni),
& aium(ni),aqum(ni),auum(ni),avum(ni),pu(ni),
& aidm(ni),aqdm(ni),audm(ni),avdm(ni),pd(ni),
& aa(9),bb(3),
& am(1),bm(1),zal(nqmax),val(nqmax),zbe(nqmax),vbe(nqmax),
& sbr(nq),sbl(nq),uu(nq),cc(nq),szero(nqmax),
& a(nq),b(nq),s(nq),
& az(4*nqmax),bz(4*nqmax),betaz(4*nqmax),ala(4),
& aul0(ni),aur0(ni),adl0(ni),adr0(ni)
stol=1.d0-1.d-14
g=(1.d0-qs1)/(1.d0+2.d0*qs1)
m1=alamb.ge.stol.and.qs1.ge.stol
call zeropol(sbr,sbl,alamb,qs1,uu,cc,nq)
pi=3.141592653589793d0
tol=1d-12
do 10 liw=1,4
miw=liw+5
if (liw.eq.1.or.liw.eq.2) ala(liw)=alamb*qs1
if (liw.eq.3.or.liw.eq.4) ala(liw)=alamb*(1.d0-3.d0*g)/
& (1.d0+2.d0*g)
alam=ala(liw)
call zeronat(szero,alam,uu,cc,nq,miw)
do 10 i=1,nq
kliw=liw+4*(i-1)
betaz(kliw)=szero(i)
10 continue
do 90 liw=1,4
miw=liw+5
do 20 i=1,nq
kliw=liw+4*(i-1)
20 s(i)=betaz(kliw)
if (tau0.lt.1.d5) call ficoeff(a,b,s,uu,alam,tau0,nq)
if (tau0.gt.1.d5) then
call sicoeff(a,s,uu,nq)
do 1122 i=1,nq
b(i)=0d0
1122 continue
endif
do 92 i=1,nq
kliw=liw+4*(i-1)
az(kliw)=a(i)
92 bz(kliw)=b(i)
90 continue
do 30 liw=1,4
miw=liw+5
alam=ala(liw)
do 40 i=1,nq
kliw=liw+4*(i-1)
s(i)=betaz(kliw)
a(i)=az(kliw)
40 b(i)=bz(kliw)
if (tau0.le.1.d5) then
call fintst(au,ad,alam,tau,tau0,uinc,us,s,a,b,ni,nq)
else
call sintst(au,ad,tau,alamb,uinc,us,s,a,ni,nq)
endif
do 100 i=1,ni
kliw=liw+4*(i-1)
auu(kliw)=au(i)*4.d0
100 add(kliw)=ad(i)*4.d0
30 continue
c
c liw=1 corresponds to fii super 1
c liw=2 corresponds to fii super 2
c liw=3 corresponds to v super 0
c liw=4 corresponds to v super 1
c
kp=3
kp2=kp*kp
c
call albega(zal,val,zbe,vbe,alamb,uinc,qs1,fl,fr,
& aa,bb,kp2,kp,sbr,sbl,uu,cc,nq)
c
if (m1) then
call coefcon(bm,tau0,uinc,alb,fl,fr,bb,kp,am,
& sbr,sbl,uu,cc,nq)
else
call coefncon(bm,alamb,tau0,uinc,alb,fl,fr,bb,kp,am,
& zal,val,zbe,vbe,sbr,sbl,uu,cc,nq)
endif
c
80 dummy=0d0
c
c
call intarblr(aiul,aidl,aiur,aidr,alamb,tau0,qs1,
& ni,us,uinc,alb,fl,fr,tau,bb,kp,bm,zal,val,zbe,vbe,
& sbr,sbl,uu,cc,nq)
c
do 110 i=1,ni
aul0(i)=aiul(i)
aur0(i)=aiur(i)
adl0(i)=aidl(i)
110 adr0(i)=aidr(i)
ui2=uinc*uinc
sinc=dsqrt(1.d0-ui2)
s1=dsin(fi)
s2=dsin(2.d0*fi)
c1=dcos(fi)
c2=dcos(2.d0*fi)
do 60 i=1,ni
v=us(i)
v2=v*v
sv=dsqrt(1.d0-v2)
vu=v*uinc
c
kliw1=1+4*(i-1)
kliw2=2+4*(i-1)
kliw3=3+4*(i-1)
kliw4=4+4*(i-1)
aul=aiul(i)+sinc*sv*.375d0*auu(kliw1)*
& (-2.d0*vu*c1*fl+v*s1*fu)+
& .1875d0*auu(kliw2)*(vu*vu*c2*fl-
& v2*c2*fr-v*vu*s2*fu)
c
adl=aidl(i)+sinc*sv*.375d0*add(kliw1)*
& (2.d0*vu*c1*fl-v*s1*fu)+
& .1875d0*add(kliw2)*(vu*vu*c2*fl-
& v2*c2*fr-v*vu*s2*fu)
c
aur=aiur(i)+.1875d0*auu(kliw2)*
& (-ui2*c2*fl+c2*fr+uinc*s2*fu)
c
adr=aidr(i)+.1875d0*add(kliw2)*
& (-ui2*c2*fl+c2*fr+uinc*s2*fu)
c
aium(i)=aul+aur
aidm(i)=adl+adr
aqum(i)=aur-aul
aqdm(i)=adr-adl
c
auum(i)=.375d0*sinc*sv*auu(kliw1)*
& (2.d0*uinc*s1*fl+c1*fu)+
& .1875d0*auu(kliw2)*(-2.d0*v*ui2*s2*fl+
& 2.d0*v*s2*fr-2.d0*vu*c2*fu)
c
audm(i)=.375d0*sinc*sv*add(kliw1)*
& (2.d0*uinc*s1*fl+c1*fu)+
& .1875d0*add(kliw2)*(2.d0*v*ui2*s2*fl-
& 2.d0*v*s2*fr+2.d0*vu*c2*fu)
c
avum(i)=.375d0*fv*(-vu*auu(kliw3)+
& sinc*sv*auu(kliw4)*c1)
c
avdm(i)=.375d0*fv*(vu*add(kliw3)+
& sinc*sv*add(kliw4)*c1)
c
anim=aium(i)
alug=aqum(i)
pu(i)=10.d0
if (anim.ge.tol) pu(i)=dsqrt(aqum(i)*aqum(i)+
& auum(i)*auum(i)+avum(i)*avum(i))/anim
anim=aidm(i)
alug=aqdm(i)
pd(i)=10.d0
if (anim.ge.tol) pd(i)=dsqrt(aqdm(i)*aqdm(i)+
& audm(i)*audm(i)+avdm(i)*avdm(i))/anim
60 continue
return
end
c CODES FOR ISOTROPIC SCATTERING NEEDED IN SUBROUTINE STOKES
c
c
subroutine sicoeff(a,s,uu,nq)
implicit real*8(a-h,o-z)
c sicoeff determines the coefficients of the approximate
c resolvent function (Sobolev's PHI)for a semi-infinite
c medium and stores them in array a
c nq is the order of the Gaussian quadrature
parameter (nqmax=150)
dimension cji(nqmax*nqmax),a(nq),s(nq),uu(nq)
ij=0
do 100 i=1,nq
bet=s(i)
do 110 j=1,nq
ij=ij+1
cji(ij)=1.d0/(1.d0-bet*uu(j))
110 continue
a(i)=1.d0/uu(i)
100 continue
call linsys(cji,a,nq,jsc)
c
return
end
subroutine hfun(h,u,s,a,nq)
implicit real*8(a-h,o-z)
c hfun determines the H-function
c at the angular variable u (u.ge.0)
c h=h(u)
c nq is the order of the Gaussian quadrature
dimension s(nq),a(nq)
if (u.lt.0.222d-15) then
h=1d0
return
endif
ui=1.d0/u
h=0.d0
do 110 i=1,nq
h=h+a(i)/(ui+s(i))
110 continue
h=h+1.d0
return
end
subroutine higif(hi,gi,tau,u,s,a,nq,ing)
implicit real*8(a-h,o-z)
c higif determines the hi- and gi-functions
c (small h and g) at the optical depth tau
c and angular variable u (u.ge.0.)
c hi=hi(tau,u), gi=gi(tau,u)
c nq is the order of the Gaussian quadrature
c if (u.lt.infinity) then ing=1
c if (u.eq.infinity) then ing=2
c
c CAUTION!
c
c if (u.eq.infinity.and.alamb.eq.1.0) only gi
c can be determined!
c in this case hi=Hopf's function
dimension s(nq),a(nq)
data tol/0.222d-15/
if (u.le.tol) then
hi=1.d0
gi=0.0
return
endif
if (ing.eq.2) then
ii=1
if (s(1).lt.tol ) ii=2
s1=0d0
s2=0d0
do 120 i=ii,nq
ai=a(i)
be=s(i)
ex=dexp(-be*tau)
ab=ai/be
s1=s1+ab
s2=s2+ab*ex
120 continue
if (s(1).gt.tol) then
hi=1.d0+s2
gi=1.d0+s1-s2
return
else
a1=a(1)
gr=1.d0+s1-s2
gi=a1*tau+gr
hi=gr/a1
return
endif
else
ui=1.d0/u
uip=ui+tol
uim=ui-tol
hi=0d0
gi=0d0
et=dexp(-tau*ui)
do 130 i=1,nq
be=s(i)
ex=dexp(-be*tau)
aa=a(i)
hi=hi+aa*ex/(ui+be)
if (be.ge.uim.and.be.le.uip) ee=aa*tau*ex
if (be.lt.uim.or.be.gt.uip) ee=aa*(ex-et)/(ui-be)
gi=gi+ee
130 continue
hi=hi+1.d0
gi=gi+et
endif
return
end
subroutine dgkap(dhu,dgu,u,tau,s,a,nq)
implicit real*8(a-h,o-z)
c dgkap computes the derivatives of
c hi- and gi-functions with respect to angular
c variable u at optical depth tau
c dhu=dhi(tau,u)/du
c dgu=dgi(tau,u)/du
dimension s(nq),a(nq)
s1=0.d0
s2=0.d0
s3=0.d0
s4=0.d0
ex=dexp(-tau/u)
do 100 i=1,nq
ai=a(i)
be=s(i)
tb=1.d0-be*u
tg=1.d0+be*u
tb2=tb*tb
exc=dexp(-be*tau)
s1=s1+ai*exc/tb2
s2=s2+ai/tb2
s3=s3+ai/tb
s4=s4+ai*exc/(tg*tg)
100 continue
dhu=s4
dgu=ex*(tau/(u*u)-s2-tau*s3/u)+s1
return
end
subroutine sintst(au,ad,tau,alamb,u,us,s,a,ni,nq)
implicit real*8(a-h,o-z)
c sintst computes the intensities of the upward (au) and
c downward (ad) moving radiation at an optical depth tau
c and at angles defined by an array us in a semi-infinite
c atmosphere (the elements of the array
c are cosines of these angles!). The parallel beam of incident
c radiation illuminates the upper boundary of the atmosphere
c at an angle arc cos (u).
dimension au(ni),ad(ni),us(ni),s(nq),a(nq)
al4=.25d0*alamb
call hfun(h,u,s,a,nq)
call higif(hu,gu,tau,u,s,a,nq,1)
do 100 iv=1,ni
v=us(iv)
call higif(hv,gv,tau,v,s,a,nq,1)
au(iv)=al4*u*h*(gu+hv-1.d0)/(u+v)
if (u.ne.v) then
ad(iv)=al4*u*h*(gv-gu)/(v-u)
else
call dgkap(dhu,dgu,u,tau,s,a,nq)
ad(iv)=al4*u*h*dgu
endif
100 continue
return
end
subroutine ficoeff(a,b,s,uu,alamb,tau0,nq)
implicit real*8(a-h,o-z)
parameter (nqmax=150)
c ficoeff determines the coefficients for an approximate resolvent
c function (SOBOLEV's PHI) in the case of an optically
c finite medium and stores them in arrays a and b
c nq is the order of the Gaussian quadrature
c
dimension cji(nqmax*nqmax),dji(nqmax*nqmax),ppc(nqmax),ppd(nqmax),
& s(nq),a(nq),b(nq),uu(nq)
data tol/.222d-15/
ii=1
if (s(1).le.tol) ii=2
do 100 j=1,nq
uj=uu(j)
uju=1.d0/uj
do 110 i=ii,nq
ij=(i-1)*nq+j
bet=s(i)
ex=dexp(-tau0*bet)
z=bet*uj
ze=(1.d0-z)*ex
z2=1.d0-z*z
cji(ij)=(1.d0+z-ze)/z2
dji(ij)=(1.d0+z+ze)/z2
110 continue
ppc(j)=uju
ppd(j)=uju
100 continue
if (ii.ne.1) then
do 120 i=1,nq
dji(i)=1.d0
cji(i)=-tau0-2.d0*uu(i)
120 continue
endif
call linsys(dji,ppd,nq,jsd)
call linsys(cji,ppc,nq,jsc)
do 140 i=ii,nq
pc=ppc(i)
pd=ppd(i)
a(i)=0.5d0*(pd+pc)
b(i)=0.5d0*(pd-pc)
140 continue
if (ii.ne.1) then
pc=ppc(1)
a(1)=0.5d0*(ppd(1)-pc*tau0)
b(1)=pc
endif
return
end
subroutine xyfun(x,y,u,tau0,s,a,b,nq,ing)
implicit real*8(a-h,o-z)
c xyfun computes the x- and y-functions (small x and small y)
c at the angular variable u
c x=x(u,tau0), y=y(u,tau0)
c if (u.lt.infinity) then ing=1
c if (u.eq.infinity) then ing=2
c nq is the order of the Gaussian quadrature
dimension s(nq),a(nq),b(nq)
data tol/0.222d-15/
ii=1
if (s(1).le.tol) ii=2
if (ing.ne.1) then
s1=0.d0
if (ii.eq.2) s1=(a(1)+.5d0*b(1)*tau0)*tau0
130 s2=0.d0
do 140 i=ii,nq
be=s(i)
ex=dexp(-be*tau0)
s2=s2+(a(i)+b(i))*(1.d0-ex)/be
140 continue
x=1.d0+s1+s2
y=x
return
endif
exy=dexp(-tau0/u)
s1=0.d0
s2=0.d0
s3=0.d0
s4=0.d0
uin=1.d0/u
uip=uin+tol
uim=uin-tol
do 110 i=ii,nq
bet=s(i)
aa=a(i)
bb=b(i)
ex=dexp(-tau0*bet)
ext=dexp(-tau0*(bet+uin))
bmz=uin-bet
bpz=uin+bet
qs=(1.d0-ext)/bpz
s1=s1+aa*qs
s4=s4+bb*qs
if (bet.ge.uim.and.bet.le.uip) then
gx=tau0*ex
s2=s2+bb*gx
s3=s3+aa*gx
else
qr=(ex-exy)/bmz
s2=s2+bb*qr
s3=s3+aa*qr
endif
110 continue
s5=0.d0
s6=0.d0
s7=0.d0
if (ii.ne.1) then
a1=a(1)
b1=b(1)
au1=a1*u
bu1=b1*u
s8=1.d0-exy
s5=au1*s8
s6=bu1*(u-exy*(tau0+u))
s7=bu1*(tau0-u*s8)
endif
x=1.d0+s5+s6+s1+s2
y=exy+s5+s7+s3+s4
return
end
subroutine xiyif(xi,yi,u,tau,tau0,s,a,b,nq,ing)
implicit real*8(a-h,o-z)
c xiyif computes the xi- and yi-functions
c at the optical depth tau and angular
c variable u
c xi=xi(tau,u,tau0), yi=yi(tau,u,tau0)
c if (u.lt.infinity) then ing=1
c if (u.eq.infinity) then ing=2
c nq is the order of the Gaussian quadrature
dimension s(nq),a(nq),b(nq)
data tol/0.954d-13/
ii=1
ta=tau0-tau
if (s(1).lt.tol) then
ii=2
a1=a(1)
b1=b(1)
endif
110 if (ing.eq.2) then
ss=0d0
si=0d0
su=0d0
sl=0d0
sd=0d0
sf=0d0
do 130 i=ii,nq
bet=s(i)
aa=a(i)
bb=b(i)
ex1=dexp(-bet*tau)
ex2=dexp(-bet*ta)
ex3=dexp(-bet*tau0)
aab=aa/bet
bbb=bb/bet
ss=ss+aab*ex1
si=si+aab*ex3
su=su+bbb*ex3
sl=sl+bbb
sd=sd+aab
sf=sf+bbb*ex2
130 continue
c1=0d0
c2=0d0
c3=0d0
c4=0d0
if (s(1).lt.tol) then
c1=a1*ta
c2=0.5d0*b1*ta*(tau+tau0)
c3=a1*tau
c4=0.5d0*b1*tau*tau
endif
xi=1.d0+c1+c2+ss-si-sf+sl
yi=1.d0+c3+c4+sf-su-ss+sd
return
endif
if (u.lt.tol) then
xi=1.d0
yi=0.0
return
endif
vu=1.d0/u
vup=vu+tol
vum=vu-tol
et=dexp(-tau*vu)
eta=dexp(-ta*vu)
c1=0d0
c2=0d0
c3=0d0
c4=0d0
do 150 i=ii,nq
bet=s(i)
aa=a(i)
bb=b(i)
ex1=dexp(-bet*tau)
ex2=dexp(-bet*tau0-vu*ta)
ex3=dexp(-bet*ta)
ex4=dexp(-bet*tau0-vu*tau)
bm=vu-bet
bp=vu+bet
aap=aa/bp
bbp=bb/bp
c1=c1+aap*(ex1-ex2)
c4=c4+bbp*(ex3-ex4)
if (bet.ge.vum.and.bet.le.vup) then
c2=c2+bb*ta*ex3
c3=c3+aa*tau*ex1
else
c2=c2+bb*(ex3-eta)/bm
c3=c3+aa*(ex1-et)/bm
endif
150 continue
c9=0d0
c10=0d0
c11=0d0
c12=0d0
if (s(1).lt.tol) then
c9=a1*u*(1.d0-eta)
c10=b1*u*(tau+u-(tau0+u)*eta)
c11=a1*u*(1.d0-et)
c12=b1*u*(tau-u+u*et)
endif
xi=1.d0+c9+c10+c1+c2
yi=et+c11+c12+c3+c4
return
end
subroutine dxidyi(dxu,dyu,tau0,tau,u,s,a,b,nq)
implicit real*8(a-h,o-z)
c dxidyi computes the derivatives of x- and
c y-functions with respect to u at optical
c depth tau
c dxu=dxi(tau,u,tau0)/du
c dyu=dyi(tau,u,tau0)/du
c nq is the order of the Gaussian quadrature
dimension s(nq),a(nq),b(nq)
data tol/.95d-13/
ii=1
if (s(1).lt.tol) ii=2
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
ta=tau0-tau
ex1=dexp(-ta/u)
ex2=dexp(-tau/u)
do 100 i=ii,nq
be=s(i)
aa=a(i)
bb=b(i)
e1=dexp(-be*ta)
e2=dexp(-be*tau)
e3=dexp(-be*tau0)
tp=1.d0+be*u
tm=1.d0-be*u
tp2=tp*tp
tm2=tm*tm
ap2=aa/tp2
bm2=bb/tm2
am2=aa/tm2
bp2=bb/tp2
s1=s1+ap2*e2
s2=s2+bm2*e1
s3=s3+am2*e2
s4=s4+bp2*e1
s5=s5+aa*e3/tp
s6=s6+bb/tm
s7=s7+aa/tm
s8=s8+bb*e3/tp
s9=s9+ap2*e3
s10=s10+bm2
s11=s11+am2
s12=s12+bp2*e3
100 continue
xa=0.d0
ya=0.d0
xb=0.d0
yb=0.d0
xc=0.d0
yc=0.d0
if (s(1).lt.tol) then
a1=a(1)
b1=b(1)
xa=a1+b1*(tau0+u)
xb=a1+b1*(tau0+2.d0*u)
ya=-a1*u+b1*u*u
yb=-a1+2.d0*b1*u
xc=a1+b1*(tau+2.d0*u)
yc=a1+b1*(tau-2.d0*u)
endif
dxu=xc+s1+s2-(ta/u*(xa+s5+s6)+
& xb+s9+s10)*ex1
dyu=yc+s3+s4+((1.d0+ya-u*(s7+s8))*tau/(u*u)+
& yb-s11-s12)*ex2
return
end
subroutine fintst(au,ad,alamb,tau,tau0,u,us,s,a,b,ni,nq)
implicit real*8(a-h,o-z)
c fintst computes the intensities of the upward (au) and
c downward (ad) moving radiation at an optical depth tau
c and at angles defined by an array us while
c the elements of the array us are cosines of these angles.
c The parallel beam of incident
c radiation illuminates the upper boundary of the atmosphere
c at an angle arc cos (u).
c The atmosphere is optically finite and its optical
c thickness is tau0.
c nq is the order of the Gaussian quadrature
dimension au(ni),ad(ni),us(ni),s(nq),a(nq),b(nq)
al4=.25d0*alamb
ta=tau0-tau
call xyfun(x,y,u,tau0,s,a,b,nq,1)
call xiyif(xu, yu, u,tau,tau0,s,a,b,nq,1)
call xiyif(xtu,ytu,u,ta, tau0,s,a,b,nq,1)
do 100 iv=1,ni
v=us(iv)
call xiyif(xv, yv, v,tau,tau0,s,a,b,nq,1)
call xiyif(xtv,ytv,v,ta, tau0,s,a,b,nq,1)
au(iv)=al4*u*(x*(xv+yu-1.d0)-
& y*(xtu+ytv-1.d0))/(u+v)
if (u.ne.v) then
ad(iv)=al4*u*(x*(yv-yu)-y*(xtv-xtu))/(v-u)
else
call dxidyi(dxv, dyv, tau0,tau,v,s,a,b,nq)
call dxidyi(dxtv,dytv,tau0,ta, v,s,a,b,nq)
ad(iv)=al4*v*(x*dyv-y*dxtv)
endif
100 continue
return
end