c      HOMOGEN.FOR
c
c A PACKAGE OF CODES TO FIND THE MONOCHROMATIC RADIATION FIELD IN
c A HOMOGENEOUS PLANE-PARALLEL ATMOSPHERE WHICH SCATTERS ACCORDING TO SOME SIMPLE PHASE FUNCTIONS
c
c
c   10.01.2010  (dd.mm.yyy)
c
c

c T.Viik
c Tartu Observatory
c 61602 Tartu maakond
c ESTONIA
c phone: 372 7 410 154
c cell phone: 372 50 89 045
c fax: 372 7 410 205
c e-mail: viik@aai.ee
c homepage: http://www.aai.ee/~viik

c  A   General information

c This package allows to find the monochromatic radiation field in a homogeneous scattering
c plane-parallel atmosphere by approximating the Sobolev's resolvent function.
c Some simple types of phase functions have been considered: isotropic scattering,
c linear anisotropic scattering, scalar Rayleigh scattering.
c
c The well-known standard problems like the planetary problem (parallel beam is incident on
c the upper boundary of an atmosphere), the Milne problem (sources of radiation are
c infinitely deep in the atmosphere), problems with different internal sources have been solved
c in an explicit way.

c The codes are written in FORTRAN-77 in double precision.

c You can find the package here - http://www.aai.ee/~viik/homogen.for

c The input parameters are:

c    - the optical thickness of the atmosphere t_0 = tau_0,

c    - optical depth in the atmosphere t = tau ,

c    - angle between the direction of photon's flight and the downward normal to
c      the upper surface of the atmosphere µ = arc cos u ,

c    - the albedo of the single scattering (or the probability of photon survival
c      in the act of interaction) \lambda = alamb ,

c    - iw  characterizes the type of scattering (isotropic, linear anisotropic,
c      scalar Rayleigh, different Mullikin indeces)

c    - xch is the parameter of linear anistropic scattering (Chandrasekhar's x)

c -    the order of Gaussian quadrature N = nq in the interval (0,1),

c -    uu ={ uu (1),..., uu ( nq )} - array of the Gaussian points, while uu (1)< uu (2)<...< uu ( nq ),
c -    cc ={ cc (1),..., cc ( nq )} - array of the respective weights,

c -    s ={ s (1),..., s ( nq )} - array of zeros of the approximate
c      characteristic equation, while s (1)< s (2)<...< s ( nq ),

c -    a ={ a (1),..., a ( nq )} - array of coefficients a i of the
c      approximate resolvent function,

c -    b ={ b (1),..., b ( nq )} - array of coefficients b i of the
c      approximate resolvent function.



c  B    LIST OF CODES

c         GENERAL CODES

c     1. Subroutine gauleg ( x 1, x 2, uu , cc , nq ) - Rybicki's subroutine from
c        ''Numerical Recipes'' to find the Gaussian points and weights in the
c        interval ( x 1, x 2);

c     2. Subroutine linsys ( a , b , n , ks ) - solves a set of linear algebraic
c        equations, a - one-dimensional array of coefficients, b - vector of
c        right-hand sides, after completion it contains the solution, n - number of
c        equations, ks - index of the result;

c     3. Subroutine zero ( s , alamb , iw, xch, uu , cc , nq ) - finds the array
c        of zeros s of the approximate characteristic equation;

c     5. Function fgx ( x , alamb , iw, xch, uu , cc , nq ) - finds the value of
c        the approximate characteristic function;

c     6. Function zeroin ( ax , bx , f , tol , alamb , iw, xch, uu , cc , nq ) -
c        zeroin is the zero of function f if found in the interval ( ax , bx ).
c        This code is for solving the APPROXIMATE characteristic equation;

c     7. Function fgak ( x , alamb ) - finds the value of the exact characteristic
c        function;

c     8. Function zeroex ( ax , bx , f , tol , alamb ) - zeroex is the zero of function f
c        if found in the interval ( ax , bx ). This code is for solving the EXACT
c        characteristic equation;


c      CODES FOR OPTICALLY SEMI-INFINITE ATMOSPHERES


c     9. Subroutine sicoeff ( a , s , uu , nq )- finds the array of coefficients a for
c        the approximate resolvent function;

c     10. Subroutine siresolv ( res , tau , taupr , s , a , nq ) - finds the resolvent res at optical depths tau and taupr ;

c     11. Subroutine hfun ( h , u , s , a , nq ) - finds the H - function h at an
c         angular argument u ;

c     12. Subroutine higif ( hi , gi , tau , u , s , a , nq , ing ) - finds the h - and g -
c         functions hi and gi at arguments tau and u ;

c     13. Subroutine dgkap ( dhu , dgu , u , tau , s , a , nq ) - finds the derivative of h -
c         and g - functions ( dhu and dgu , respectively) with respect to u at
c         arguments tau and u ;

c     14. d2gkap -  d2gkap computes the second derivative of gi-functions with respect to angular
c         variable u at optical depth tau -  d2gu=d2gi(tau,u)/du2

c     15. Subroutine sintst ( au , ad , tau , alamb , u , us , s , a , ni , nq ) - for planetary
c         problem finds arrays of the intensities of upward and downward moving
c         radiation au and ad at the incident angle arc cos u and at the optical
c         depth tau . The intensities are found at given angles in array us ={ us (1),..., us ( ni )};

c     16. Subroutine silin ( au , ad , alamb , tau , p 0, p 1, s , a , uu , nq ) - for problem with
c         linear internal sources B 0 ( tau )= p 0 + p 1*tau
c         silin finds the arrays of the intensities of upward and downward moving radiation au and ad at
c         optical depth tau and at Gaussian points uu .
c         ONLY FOR NONCONSERVATIVE ATMOSPHERES!

c     17. Subroutine siquadr (au, ad, sb, alamb, tau, p2, s, a, us, ni, nq) - for problems with
c         internal sources of the type B_0( tau )= p2 *tau * tau.
c         siquadr finds the arrays of the intensities of upward and downward moving radiation au and ad at
c         optical depth tau and at given points in an array us(ni) .
c         sb is the source function

c     18. Subroutine emixed(au, ad, sb, tau, alamb, u, us, s, a, ni, nq) - emixed computes the intensities
c         of the upward (au) and downward (ad) moving radiation at an optical depth tau
c         and at angles defined by an array us in a semi-infinite atmosphere (the elements of the array
c         are cosines of these angles!).
c         The internal sources are distributed according to law: B_0(tau)=tau*exp(-tau/u)
c         sb is the source function

c     19. Subroutine blayer(au, ad, sb, tau, b0, tau1, s, a, us, ni, nq)
c         Blayer finds the intensities au and ad and the source function sb in a homogeneous
c         isotropic semi-infinite atmosphere at angles defined by an array us(ni) in a
c         semi-infinite atmosphere (the elements of the array
c         are cosines of these angles!) where there is a constant source of strength b0
c         in a infinitesimally thin layer at tau=tau_1 : B_0 (tau)= b0 * delta( tau - tau_1 ).


c     20. Subroutine milne(au, ad, sb, tau, alamb, s, a, us, ni, nq)
c         Milne computes the intensities in a homogeneous  scattering semi-infinite medium
c         (with sources at infinity) at optical depth tau and at arbitrary angular points, defined by an
c         array us(ni).
c         sb is the source function normalized to sb ( tau = 0 )=1 .



c      CODES FOR OPTICALLY FINITE ATMOSPHERESc



c     21. Subroutine ficoeff ( a , b , s , uu , alamb , tau 0, nq ) - finds the arrays of
c         coefficients a and b for the approximate resolvent function

c     22. Subroutine firesolv ( res , tau 0, tau , taupr , s , a , b , nq ) - finds the resolvent
c         res at optical depths tau and taupr

c     23. Subroutine xyfun ( x , y , u , tau 0, s , a , b , nq , ing ) - finds the X- and Y-functions
c         at the angular variable u.
c         If (u.lt.infinity) then ing=1
c         if (u.eq.infinity) then ing=2


c     24. Subroutine xiyif ( xi , yi , u , tau , tau 0, s , a , b , nq , ing ) - finds the small x - and
c         small y - functions an angle arc cos u and at optical depth tau

c     25. Subroutine dxidyi ( dxu , dyu , tau0, tau , u , s , a , b , nq ) - finds the derivative
c         of small x - and small y - functions ( dxu and dyu , respectively) with respect to u
c         at arguments tau and u

c     26. Subroutine fintst ( au , ad , alamb , tau , tau 0, u , us , s , a , b , ni , nq ) -
c         finds arrays of the intensities of upward and downward  moving radiation au and ad at the incident
c         angle arc cos u and at the optical depth tau  for a planetary problem .
c         The intensities are found at given angles in array us ={ us (1),..., us ( ni )}

c     27. Subroutine filin ( au , ad , alamb , tau 0, tau , p 0, p 1, s , a , b , uu , nq ) - filin computes the
c         intensities au and ad in a homogeneous isotropically scattering optically finite medium
c         at optical depth tau and at the points of gaussian quadrature on interval (0,1)
c         if there are internal sources B_( tau ) = p0 + p1 * tau in the medium .

c     28. Subroutine emixedf( au, ad, sb, tau, tau0, akap, us, s, a, b, ni, nq) - emixedf computes the
c         intensities au and ad of upward and downward  moving radiation in a homogeneous isotropically
c         scattering optically finite medium at optical depth tau and at the points of a given angles in array us(i),i=1,...,ni
c         if there are internal sources B_( tau ) = tau * exp(-tau/akap) in the medium .

c     29. Subroutine fisquare_out( au, ad, tau0, s, a, b, us, ni, nq) - computes the intensities au and ad of the emergent
c         radiation in a homogeneous isotropically scattering optically finite medium
c         at given angles whose cosines are in array us ={ us (1),..., us ( ni )}
c         if there are internal sources B_( tau )= tau*tau in the medium;

c     30. Subroutine xlq( q, qx5, qt5, u, si, tau0, kl) - finds the values of q, qx5 and qt5 for subroutine fisquare_out.

c     31. Subroutine qfun ( q , tau0, r , tau , s , a , b , nq ) - finds the Q - function q for subroutine filin .



c
c      3. GENERAL CODES
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 weights 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



      subroutine zero(s,alamb,iw,xch,uu,cc,nq)
      implicit real*8(a-h,o-z)
c     zero evaluates the zeros of the
c     approximate characteristic equation
c     and stores them in array s
      dimension uu(nq),s(nq),cc(nq)
      external fgx
      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=tol
      stol=1.d0-tol
      do 100 ij=1,nq
         if (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)=zeroin(ax,bx,fgx,tolex,alamb,iw,xch,uu,cc,nq)
  100    continue
c
c     zeroin = 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
c
      if (alamb.gt.stol) s(1)=0.d0
      return
      end






      function fgx(x,alamb,iw,xch,uu,cc,nq)
      implicit real*8(a-h,o-z)
c     fgx evaluates the approximate characteristic equation
c     as a function of x
c     x - argument on the beta-axis
c     alamb - probability of photon survival
c     iw - indicates the type of scattering
c     xch - parameter for linear anisotropic scattering
      dimension uu(nq),cc(nq)
      s=0.d0
      do 100 i=1,nq
      u=uu(i)
      u2=u*u
c
c     iw is the characteristic function for
c     iw=0 - isotropic scattering
c     iw=1 - linear anisotropic scattering m=0,
c     iw=2 - linear anisotropic scattering m=1,
c     iw=3 - Rayleigh scalar scattering
c     iw=4 - Rayleigh scattering m=0,
c     iw=5 - Rayleigh scattering m=1,
c     iw=6 - Rayleigh scattering m=2.
c
c     Rayleigh scattering (with Mullikin indeces)
c
c     iw=7 -  Mullikin index 1
c     iw=8 -  Mullikin index 2
c     iw=9 -  Mullikin index 3
c     iw=10 -  Mullikin index 4
c     iw=11 - Mullikin index 5
c
c
      u2=u*u
      uv=1.d0-u2
      al=1.d0-alamb
      if (iw.eq.0) w=.5d0*alamb
      if (iw.eq.1) w=.5d0*alamb*(1.d0+xch*al*u2)
      if (iw.eq.2) w=.25d0*alamb*xch*uv
      if (iw.eq.3) w=.1875d0*alamb*(3d0-u2)
      if (iw.eq.4) w=.1875d0*alamb*(3.d0*al*u2*u2-
     &     (2.d0-alamb)*u2+3.d0)
      if (iw.eq.5) w=.375d0*alamb*u2*uv
      if (iw.eq.6) w=9.375d-2*alamb*uv*uv
      if (iw.eq.7) w=0.375d0*alamb*uv*(1.d0+2.d0*u2)
      if (iw.eq.8) w=0.1875d0*alamb*(1.d0+u2)*(1.d0+u2)
      if (iw.eq.9) w=0.75d0*alamb*u2
      if (iw.eq.10) w=0.375d0*alamb*uv
      if (iw.eq.11) w=0.75d0*alamb*uv
c
      sx=1.d0/(1.d0+x*u)+1.d0/(1.d0-x*u)
      s=s+cc(i)*w*.5d0*sx
  100 continue
      fgx=1.d0-2.d0*s
      return
      end



      function zeroin(ax,bx,f,tol,alamb,iw,xch,uu,cc,nq)
      implicit real*8(a-h,o-z)
c     zeroin 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,iw,xch,uu,cc,nq)
      fb=f(b,alamb,iw,xch,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,iw,xch,uu,cc,nq)
      if ((fb*(fc/dabs(fc))).gt.0.d0) go to 20
      go to 30
   90 zeroin=b
      return
      end



      function fgak(x,alamb)
      implicit real*8(a-h,o-z)
c     this code evaluates the exact characteristic
c     equation as a function of x
      fgak=1d0-alamb*dlog((1d0+x)/(1d0-x))/(2d0*x)
      return
      end



      function zeroex(ax,bx,f,tol,om)
      implicit real*8(a-h,o-z)
c     zeroex determines zero of function f if found in the
c      interval ax,bx
c     tol - desired length of the interval
c     of uncertainty for the zero zeroex
c     G.E.Forsythe, M.A.Malcolm, C.B.Moler
c     Viik   sept   1980
      data eps/.222d-15/
      a=ax
      b=bx
      fa=f(a,om)
      fb=f(b,om)
   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,om)
      if ((fb*(fc/dabs(fc))).gt.0.d0) go to 20
      go to 30
   90 zeroex=b
      return
      end
c
c
c
c
c     4. CODES FOR OPTICALLY SEMI-INFINITE ATMOSPHERES
c
c
c
c
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=301)
      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 siresolv(res,tau,taupr,s,a,nq)
      implicit real*8(a-h,o-z)
      parameter (nqmax=301)
c     siresolv finds the resolvent function res (Sobolev's GAMMA)
c     for a semi-infinite homogeneous medium at optical
c     depths tau and taupr
      dimension etau(nqmax),eta(nqmax),etaupr(nqmax)
      dimension s(nq),a(nq)
      sq=dsqrt(3d0)
      ii=1
      if (s(1).eq.0d0) ii=2
      ta=dabs(tau-taupr)
      do 110 i=1,nq
         be=s(i)
         etau(i)=dexp(-be*tau)
         etaupr(i)=dexp(-be*taupr)
         eta(i)=dexp(-be*ta)
  110    continue
      sfi=0d0
      do 100 i=ii,nq
         sfi=sfi+a(i)*eta(i)
  100    continue
      sf=0d0
      si=0d0
      do 200 i=ii,nq
         ai=a(i)
         bei=s(i)
         sd=0d0
         su=0d0
         do 300 j=ii,nq
            aj=a(j)
            bej=s(j)
            bij=bei+bej
            sd=sd+aj*eta(j)/bij
            su=su+aj*etaupr(j)/bij
  300       continue
         sf=sf+ai*sd
         si=si+ai*etau(i)*su
  200    continue
      sx=sf-si
      if (s(1).eq.0d0) then
         sif=0d0
         do 400 i=ii,nq
            bei=s(i)
            sif=sif+a(i)*(1d0-etau(i)-
     &          etaupr(i)+eta(i))/bei
  400       continue
         te=tau
         if (tau.ge.taupr) te=taupr
         res=sq+sfi+3d0*te+sq*sif+sx
         return
      endif
      res=sfi+sx
      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 d2gkap(d2gu,u,tau,s,a,nq)
      implicit real*8(a-h,o-z)
c     d2gkap computes the second derivative of
c     gi-functions with respect to angular
c     variable u at optical depth tau
c     d2gu=d2gi(tau,u)/du2
c
      dimension s(nq),a(nq)
      s1=0.d0
      s2=0.d0
      s3=0.d0
      ex=dexp(-tau/u)
      do 100 i=1,nq
         ai=a(i)
         si=s(i)
         ts=1.d0-si*u
         ts2=ts*ts
         exc=dexp(-si*tau)
         s1=s1+ai*exc/ts2
         s2=s2+ai/ts2
         s3=s3+ai/ts
  100    continue
      su=tau/(u*u)-s2-tau*s3/u
      z1=tau/(u*u)*ex*su
      z2=0d0
      z4=0d0
      z5=0d0
      do 110 i=1,nq
         ai=a(i)
         si=s(i)
         ts=1d0-si*u
         ts2=ts*ts
         ts3=ts2*ts
         exc=dexp(-si*tau)
         z2=z2+ai*si/ts3
         z4=z4+ai*si/ts2
         z5=z5+ai*si*exc/ts3
  110    continue
      d2gu=z1+ex*(-2d0*tau/(u*u*u)-2d0*z2+tau/(u*u)*s3-tau/u*z4)
     &       +2d0*z5
      return
      end



      subroutine sintst(au,ad,sb,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)
      sb=al4*hu*gu
      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.eq.v) go to 110
         ad(iv)=al4*u*h*(gv-gu)/(v-u)
         go to 100
  110    call dgkap(dhu,dgu,u,tau,s,a,nq)
         ad(iv)=al4*u*h*dgu
  100 continue
      return
      end



      subroutine silin(au,ad,sb,alamb,tau,p0,p1,s,a,us,ni,nq)
      implicit real*8(a-h,o-z)
c     silin computes the intensities in a homogeneous
c     isotropically scattering semi-infinite medium
c     at optical depth tau and at arbitrary points given by
c     array us(ni) if there are internal sources
c     b_0(tau)=p0+p1*tau in the medium
c     au(i) - upward intensity
c     ad(i) - downward intensity
C     sb is the source function
c
c     caution!!! this code is to be used only if lambda.lt.1.0
c
c
      dimension au(ni),ad(ni),s(nq),a(nq),uu(nq),us(ni)
      hinf=1d0/dsqrt(1d0-alamb)
      call higif(hf,ginf,tau,1d0,s,a,nq,2)
      sb0=p0*hinf*ginf
      dh=0d0
      dhi=0d0
      do 100 i=1,nq
         ai=a(i)
         si=s(i)
         si2=si*si
         dh=dh+ai/si2
         dhi=dhi+ai*dexp(-si*tau)/si2
  100    continue
      sb1=p1*(tau*hinf+ginf-dh+dhi)
      sb=sb0+sb1
      dg=hinf*tau-dh+dhi
      do 110 i=1,ni
         ui=us(i)
         call higif(hu,gu,tau,ui,s,a,nq,1)
         ar=ginf+hu-1d0
         at=ginf-gu
         au(i)=p0*hinf*ar+p1*((dh+ui*hinf)*ar+hinf*dg)
         ad(i)=p0*hinf*at+p1*((dh-ui*hinf)*at+hinf*dg)
  110    continue
      return
      end



      subroutine siquadr(au,ad,sb,alamb,tau,p2,s,a,us,ni,nq)
      implicit real*8(a-h,o-z)
c     siquadq computes the intensities in a homogeneous
c     isotropically scattering semi-infinite medium
c     at optical depth tau and at the points of an array us(ni)
c     if there is internal source
c     b_0(tau)=p2*tau*tau in the medium
c     au(i) - upward intensity
c     ad(i) - downward intensity
c     sb - source function
c
c     CAUTION!!! this code is to be used only if lambda.lt.1.0
c
c
      dimension au(ni),ad(ni),s(nq),a(nq),us(ni)
      alx=1d0/(1d0-alamb)
      hinf=1d0/dsqrt(1d0-alamb)
      sr=0d0
      sp=0d0
      s1=0d0
      s2=0d0
      do 100 i=1,nq
         ai=a(i)
         si=s(i)
         si2=si*si
         si3=si2*si
         ex=dexp(-si*tau)
         sr=sr+ai/si2
         sp=sp+ai/si3
         s1=s1+ai/si2*ex
         s2=s2+ai/si3*ex
  100    continue
      call higif(hf,ginf,tau,1d0,s,a,nq,2)
      sb=2d0*p2*(sp*ginf+tau*tau/2d0*alx-sr*sr+sp*hinf+sr*s1-hinf*s2)
c
      do 110 i=1,ni
         ui=us(i)
         call higif(hu,gu,tau,ui,s,a,nq,1)
         s3=0d0
         s4=0d0
         s5=0d0
         s6=0d0
         exu=dexp(-tau/ui)
         exu1=1d0-exu
         do 120 k=1,nq
            ak=a(k)
            sk=s(k)
            sk2=sk*sk
            sk3=sk2*sk
            ex=dexp(-sk*tau)
            s3=s3+ak/sk2*ex/(1d0+sk*ui)
            s4=s4+ak/sk3*ex/(1d0+sk*ui)
            if (sk*ui.ne.1d0) then
               s5=s5+ak/sk2*(ex-exu)/(1d0-sk*ui)
               s6=s6+ak/sk3*(ex-exu)/(1d0-sk*ui)
            else
               s5=s5+ak/sk2*tau*ex/ui
               s6=s6+ak/sk3*tau*ex/ui
            endif
  120       continue
         up=2d0*sp*(ginf+hu-1d0)-2d0*sr*sr+2d0*hinf*sp
         up=up+2d0*sr*s3-2d0*hinf*s4
         up=up+alx*(tau*tau+2d0*tau*ui+2d0*ui*ui)
         au(i)=p2*up
         down=2d0*sp*(ginf-gu)-2d0*sr*sr*exu1
         down=down+2d0*sr*s5-2d0*hinf*s6+2d0*hinf*sp*exu1
         down=down+alx*(tau*tau-2d0*ui*tau+2d0*ui*ui*exu1)
         ad(i)=p2*down
  110    continue
      return
      end



      subroutine emixed(au,ad,sb,tau,u,us,s,a,ni,nq)
      implicit real*8(a-h,o-z)
c     emixed 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!).
c     The internal sources are distributed according to law:
c     B_0(t)=t*exp(-t/u)
c     sb is the source function
c
c      Viik  03 03 2007
c
      dimension au(ni),ad(ni),us(ni),s(nq),a(nq)
      u2=u*u
      s1=0d0
      do 10 i=1,nq
         ai=a(i)
         si=s(i)
         ts=1d0+si*u
         ts2=ts*ts
         s1=s1+ai/ts2
   10    continue
      dh=s1
      call hfun(h,u,s,a,nq)
      call higif(hu,gu,tau,u,s,a,nq,1)
      call dgkap(dhu,dgu,u,tau,s,a,nq)
      sb=u2*(dgu*h+gu*dh)
c
      do 20 i=1,ni
         ui=us(i)
         call higif(hi,gi,tau,ui,s,a,nq,1)
         z=gu+hi-1d0
         au(i)=u2*(h/(u+ui)*z+u/(u+ui)*dh*z-u*h/(u+ui)/(u+ui)*z+
     &         u*h/(u+ui)*dgu)
         if (u.ne.ui) then
            ad(i)=u2*(h*(gu-gi)/(u-ui)+u/(u-ui)*dh*(gu-gi)-
     &            u*h/(u-ui)/(u-ui)*(gu-gi)+u*h/(u-ui)*dgu)
         else
            call d2gkap(d2gu,u,tau,s,a,nq)
            ad(i)=u2*(h*dgu+u*dh*dgu+.5d0*u*h*d2gu)
         endif
   20    continue
c
      return
      end



         subroutine blayer(au,ad,sb,tau,b0,tau1,s,a,us,ni,nq)
         implicit real*8(a-h,o-z)
c
c        this code finds the intensities au and ad and
c        the source function sb in a homogeneous
c        isotropic semi-infinite atmosphere where there is a
c        constant source of strength b0 in a infinitesimally
c        thin layer at tau=tau_1.
c
c        CAUTION!!! The atmosphere is non-conservative!
c
c        Viik    18 04 2007
c
         dimension a(nq),s(nq),us(ni),au(ni),ad(ni)
         do 100 k=1,ni
            u=us(k)
            call higif(hiu,giu,tau,u,s,a,nq,1)
            call higif(hiu1,giu1,tau1,u,s,a,nq,1)
            if (tau.le.tau1) then
               dxt=dexp(-(tau1-tau)/u)
               s1=0d0
               s3=0d0
               s5=0d0
               s7=0d0
               s8=0d0
               s10=0d0
               s11=0d0
               do 110 i=1,nq
                  ai=a(i)
                  si=s(i)
                  ps=1d0/si
                  dxs=dexp(-si*(tau1-tau))
                  call hfun(h,ps,s,a,nq)
                  call higif(hi,gi,tau,ps,s,a,nq,1)
                  call higif(hi1,gi1,tau1,ps,s,a,nq,1)
                  denm=1d0/(1d0-si*u)
                  denp=1d0/(1d0+si*u)
                  s1=s1+ai*h*denm*(dxs-dxt)
                  s3=s3+ai*dexp(-si*tau1)*denm*(hi-hiu)
                  s5=s5+ai*dexp(-si*tau1)*denm*(hi1-hiu1)
                  s7=s7+ai*h*denp
                  s8=s8+ai*dexp(-si*tau1)*denp*(hi1-1d0)
                  s10=s10+ai*h*dxs*denp
                  s11=s11+ai*dexp(-si*tau1)*denp*
     &                   (hi+giu-1d0)
  110             continue
            au(k)=b0*(s1-s3+dxt*(s5+s7-s8)+dxt/u)
            ad(k)=b0*(s10-s11)
         else
               dxt=dexp(-(tau-tau1)/u)
               s12=0d0
               s13=0d0
               s14=0d0
               s15=0d0
               s16=0d0
               s17=0d0
               do 120 i=1,nq
                  ai=a(i)
                  si=s(i)
                  ps=1d0/si
                  call hfun(h,ps,s,a,nq)
                  call higif(hi,gi,tau,ps,s,a,nq,1)
                  call higif(hi1,gi1,tau1,ps,s,a,nq,1)
                  denm=1d0/(1d0-si*u)
                  denp=1d0/(1d0+si*u)
                  s12=s12+ai*h*denp
                  s13=s13+ai*dexp(-si*tau1)*denp*(hi1+giu1-1d0)
                  s14=s14+ai*h*denm*(dexp(-si*(tau-tau1))-dxt)
                  s15=s15+ai*denm*(hi1-1d0)*
     &                (dexp(-si*tau)-dexp(-si*tau1)*dxt)
                  s16=s16+ai*h*denp*dexp(-si*(tau-tau1))
                  s17=s17+ai*dexp(-si*tau)*denp*(hi1-1d0)
  120             continue
            au(k)=b0*(s16-s17)
            ad(k)=b0*(dxt*s12-dxt*s13+s14-s15+dxt/u)
         endif
  100    continue
      call siresolv(sb,tau,tau1,s,a,nq)
      sb=b0*sb
      return
      end



      subroutine milne(au,ad,sb,tau,alamb,s,a,us,ni,nq)
      implicit real*8(a-h,o-z)
c     milne computes the intensities in a homogeneous
c     isotropically scattering semi-infinite medium
c     (with sources at infinity) at optical depth tau and
c     at arbitrary angular points, defined by an array us(ni)
c     au(i) - upward intensity
c     ad(i) - downward intensity
      dimension au(ni),ad(ni),s(nq),a(nq),us(ni)
      external fgak
      data tol/.954d-15/
      if (s(1).lt.tol) then
         call higif(q,gi,tau,1d0,s,a,nq,2)
         sqt=(tau+q)/q
         sb=sqt
         do 120 i=1,ni
            ui=us(i)
            call higif(hu,gu,tau,ui,s,a,nq,1)
            au(i)=sqt+hu-1d0
            ad(i)=sqt-gu
  120       continue
         return
      else
         ax=tol
         bx=1d0-tol
         ak=zeroex(ax,bx,fgak,tol,alamb)
         aki=1d0/ak
         call hfun(haki,aki,s,a,nq)
         call higif(hiaki,giaki,tau,aki,s,a,nq,1)
         hakex=haki*dexp(ak*tau)-hiaki
         sb=hakex+1d0
         do 130 i=1,ni
            ui=us(i)
            call higif(hu,gu,tau,ui,s,a,nq,1)
            aku=ak*ui
            au(i)=(hakex+hu)/(1d0-aku)
            ad(i)=(hakex+1d0-gu)/(1d0+aku)
  130       continue
      endif
      return
      end






c
c
c
c
c      5. CODES FOR OPTICALLY FINITE ATMOSPHERES
c
c
c
c
c
c




      subroutine ficoeff(a,b,s,uu,alamb,tau0,nq)
      implicit real*8(a-h,o-z)
      parameter (nqmax=301)
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 firesolv(res,tau0,tau,taupr,s,a,b,nq)
      implicit real*8(a-h,o-z)
      parameter (nqmax=301)
c     firesolv finds the resolvent res for an optically finite
c     homogeneous medium of optical thickness tau0
c     at optical depths tau and taupr (Sobolev's GAMMA)
c
      dimension ey(nqmax),etpt(nqmax),etot(nqmax),etotpt(nqmax),
     &          ex(nqmax),etau0(nqmax),etotp(nqmax),
     &          s(nq),a(nq),b(nq)
      data eps/.222d-15/
      ii=1
      if (s(1).lt.eps) then
         ii=2
         a1=a(1)
         b1=b(1)
      endif
      x=tau
      y=taupr
      if ((y-x).lt.0d0) then
         x=taupr
         y=tau
      endif
      s1=0d0
      tpt=y-x
      tot=tau0-x
      totp=tau0-y
      totpt=tau0-tpt
      do 150 i=1,nq
         bei=s(i)
         ey(i)=dexp(-bei*y)
         etpt(i)=dexp(-bei*tpt)
         etot(i)=dexp(-bei*tot)
         etotp(i)=dexp(-bei*totp)
         etotpt(i)=dexp(-bei*totpt)
         ex(i)=dexp(-bei*x)
         etau0(i)=dexp(-bei*tau0)
  150    continue
      if (s(1).lt.eps) then
         sq=b1*x*totp*(2d0*a1+b1*tau0)
         sf=a1+b1*x
         sd=a1+b1*tot
         su=a1+b1*y
         sl=a1+b1*totp
         s9=a1+b1*tpt
      endif
      s3=0d0
      do 140 in=ii,nq
         be=s(in)
         e1=etpt(in)
         e2=etotpt(in)
         s1=s1+a(in)*e1+b(in)*e2
  140    continue
      s4=0d0
      s5=0d0
      s6=0d0
      s7=0d0
      s8=0d0
      s10=0d0
      if (s(1).lt.eps) then
         do 160 in=ii,nq
            ai=a(in)
            bi=b(in)
            be=s(in)
            eb=1d0/be
            ec=eb*eb
            am1=etpt(in)
            am2=ey(in)
            am3=ex(in)
            am4=etotp(in)
            am5=etotpt(in)
            am6=etot(in)
            am7=etau0(in)
            h1=eb*(am1-am2)
            h2=eb*(1d0-am3)
            h3=eb*(am4-am5)
            h4=eb*(am6-am7)
            c1=be*x-1d0
            c2=be*x+1d0
            h5=ec*am1*(c1+am3)
            h6=ec*(c1+am3)
            h7=ec*am4*(1d0-am3*c2)
            h8=ec*am6*(1d0-am3*c2)
            s4=s4+(ai*sf-bi*sd)*h1
            s5=s5+(bi*sf-ai*sd)*h3
            s6=s6+(ai*su-bi*sl)*h2
            s7=s7+(bi*su-ai*sl)*h4
            s8=s8+(ai+bi)*(h5+h6+h7+h8)
  160       continue
            s10=s9-sq+s4+s5+s6+s7-b1*s8
      endif
      do 170 in=ii,nq
         bei=s(in)
         ai=a(in)
         bi=b(in)
         p3=etot(in)
         p4=ex(in)
         p5=etau0(in)
         do 180 jn=ii,nq
            bej=s(jn)
            aj=a(jn)
            bj=b(jn)
            bm=bei-bej
            bp=bei+bej
            e1=etpt(jn)
            p1=ey(jn)
            p2=etotp(jn)
            p6=etotpt(jn)
            e2=p1*p4
            e3=p2*p3
            e4=p5*p6
            e5=p5*e1
            e6=p3*p1
            e7=p2*p4
            s2=s2+(ai*aj-bi*bj)*(e1-e2-e3+e4)/bp
            if (in.ne.jn) s3=s3+(ai*bj-bi*aj)*(e5-e6-e7+p6)/bm
  180       continue
  170   continue
      res=s1+s2+s3+s10
      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
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.eq.v) go to 110
         ad(iv)=al4*u*(x*(yv-yu)-y*(xtv-xtu))/(v-u)
         go to 100
  110    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)
  100    continue
         end



      subroutine filin(au,ad,alamb,tau0,tau,p0,p1,s,a,b,uu,nq)
      implicit real*8(a-h,o-z)
c     filin computes intensities in a homogeneous
c     isotropically scattering optically finite medium
c     at optical depth tau and at the points of
c     gaussian quadrature on interval (0,1)
c     if there are internal sources b0(tau)=p0+p1*tau in the medium
c     au - upward intensity
c     ad - downward intensity
c     nq is the order of the Gaussian quadrature
      dimension au(nq),ad(nq),s(nq),a(nq),b(nq),uu(nq)
      data tol/0.954d-13/
      al1=1.d0-tol
      ta=tau0-tau
      call qfun(qt,tau0,0d0,tau,s,a,b,nq)
      call qfun(qtf,tau0,0d0,tau0,s,a,b,nq)
      call qfun(qta,tau0,ta,tau0,s,a,b,nq)
      call xiyif(xz,yf,1d0,tau,tau0,s,a,b,nq,2)
      call xiyif(xf,yz,1d0,ta, tau0,s,a,b,nq,2)
      call xyfun(xif,yif,1d0,tau0,s,a,b,nq,2)
      tfx=tau0*xif-qtf
      tx=xif*(tau*yf+ta*(xf-1d0)-qt-qta)
      do 120 i=1,nq
         ui=uu(i)
         call xiyif(xt, yt, ui,tau,tau0,s,a,b,nq,1)
         call xiyif(xta,yta,ui,ta, tau0,s,a,b,nq,1)
         ar=xif*(xt+yf-xf-yta)
         at=xif*(-yt+yf+xta-xf)
         au(i)=p0*ar+p1*(ui*ar+qtf*(xt+yf-1d0)-
     &         tfx*(xf+yta-1d0)+tx)
         ad(i)=p0*at+p1*(-ui*at-qtf*(yt-yf)+
     &         tfx*(xta-xf)+tx)
  120    continue
      return
      end


      subroutine emixedf(au,ad,sb,tau,tau0,akap,us,s,a,b,ni,nq)
      implicit real*8(a-h,o-z)
c     emixedf 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 an optically finite
c     atmosphere (the elements of the array
c     are cosines of these angles!).
c     The internal sources are distributed according to law:
c     B_0(t)=t*exp(-t/akap)
c     sb is the source function
c
c     Nota bene! Akap should not be equal to any value of the array us!!!
c
c
c
c      Viik  07 01 2010
c
      dimension au(ni),ad(ni),us(ni),s(nq),a(nq),b(nq)
c
      ing=1
      ta=tau0-tau
      call xyfun(x,y,akap,tau0,s,a,b,nq,ing)
      call xiyif(xa,ya,akap,tau,tau0,s,a,b,nq,ing)
      call xiyif(xata,yata,akap,ta,tau0,s,a,b,nq,ing)
      call dxidyi(dxa,dya,tau0,tau,akap,s,a,b,nq)
      call dxidyi(dxata,dyata,tau0,ta,akap,s,a,b,nq)
      call dxidyi(dx,dyaa,tau0,0d0,akap,s,a,b,nq)
      call dxidyi(dxaa,dy,tau0,tau0,akap,s,a,b,nq)
      ak2=akap*akap
      ak3=ak2*akap
      sb=ak2*(dx*ya+x*dya-dy*(xata-1d0)-y*dxata)
      do 10 i=1,ni
         u=us(i)
         akp=u+akap
         akm=u-akap
         akp2=akp*akp
         akm2=akm*akm
         call xiyif(xu,yu,u,tau,tau0,s,a,b,nq,ing)
         call xiyif(xuta,yuta,u,ta,tau0,s,a,b,nq,ing)
         r1=xu+ya-1d0
         r2=xata+yuta-1d0
         t=yu-ya
         w=xuta-xata
         z1=x*r1-y*r2
         z2=x*t-y*w
         au(i)=ak2/akp*(u/akp*z1+akap*(dx*r1+x*dya-dy*r2-y*dxata))
c
         ad(i)=ak2/akm*(u/akm*z2+akap*(dx*t-x*dya-dy*w+y*dxata))
c         au(i)=ak3/akp*(w1-w2)+ak2*x*(u/akp2*w3+akap/akp*dya)-
c     &         ak2*y*(u/akp2*w4-akap/akp*dxata)
c         ad(i)=ak3/akm*(w5-w6)+ak2*x*(u/akm2*w7+akap/akm*dya)-
c     &         ak2*y*(u/akm2*w8-akap/akm*dxata)
   10    continue
      return
      end


      subroutine fisquare_out(au,ad,tau0,s,a,b,us,ni,nq)
      implicit real*8(a-h,o-z)
c
c     fisquare_out computes the intensities au and ad of the emergent
c     radiation in a homogeneous isotropically scattering optically finite medium
c     at the points defined by the array us (i), i=1,...,ni
c     if there are internal sources B_( tau )= tau*tau in the medium.
c
c     au(i) - intensity of the emergent radiation at tau=0
c     ad(i) - intensity of the emergent radiation at tau=tau_0
c     us(i) - array of angles at which the intensities will be computed
c     us ={us(1),...,us(ni)}
c     nq - order of Gauss quadrature
c
c
c
      dimension au(ni),ad(ni),s(nq),a(nq),b(nq),us(ni)
      t=tau0
      t2=t*t
      t3=t*t2
      td3=t3/3d0
      t4=t3*t
      if (s(1).ge.1d-12) then
         istart=1
         a1=0d0
         b1=0d0
      else
         a1=a(1)
         b1=b(1)
         istart=2
      endif
      do 10 i=1,ni
         u=us(i)
         sxx=0d0
         call xlq(q1,qx5,qt5,u,sxx,tau0,1)
         call xlq(q4,qx5,qt5,u,sxx,tau0,4)
         s1=0d0
         s2=0d0
         s3=0d0
         s4=0d0
         s5=0d0
         s6=0d0
         s7=0d0
         s8=0d0
         do 20 k=istart,nq
            sk=s(k)
            ak=a(k)
            bk=b(k)
            um=1d0-sk*u
            up=1d0+sk*u
            dex=dexp(-sk*t)
            call xlq(q2,qx5,qt5,u,sk,tau0,2)
            call xlq(q3,qx5,qt5,u,sk,tau0,3)
            s2=s2+bk*(q3-dex*q1)/up
            s3=s3+ak*(q3-dex*q1)/up
            s6=s6+bk*(q2-dex*q4)/up
            s7=s7+ak*(q2-dex*q4)/up
            if (sk*u.ne.1d0) then
               qx1=(q2-q1)/um
               qx2=(q3-q4)/um
            else
               call xlq(q,qx5,qt5,u,sk,tau0,5)
               qx1=qx5
               qx2=qt5
            end if
               s1=s1+ak*qx1
               s4=s4+bk*qx1
               s5=s5+ak*qx2
               s8=s8+bk*qx2
   20       continue
         if (s(1).le.1d-12) then
            sum1=q1+a1*u*(td3-q1)+b1*u*(t4/4d0-u*td3+u*q1)
            sum2=a1*u*(td3-q1)+b1*u*(t4/12d0+td3*u-(t+u)*q1)
            sum3=q4+a1*u*(td3-q4)+b1*u*(t4/12d0-td3*u+u*q4)
            sum4=a1*u*(td3-q4)+b1*u*(t4/4d0+u*td3-(t+u)*q4)
         else
            sum1=q1
            sum2=0d0
            sum3=q4
            sum4=0d0
         end if
         ai1=sum1+u*(s1+s2)
         ai2=sum2+u*(s3+s4)
         ai3=sum3+u*(s5+s6)
         ai4=sum4+u*(s7+s8)
         call xyfun(x,y,u,t,s,a,b,nq,1)
         au(i)=(x*ai1-y*ai2)/u
         ad(i)=(x*ai3-y*ai4)/u
   10    continue
      return
      end


      subroutine xlq(q,qx5,qt5,u,si,tau0,kl)
      implicit real*8(a-h,o-z)
c
c     this code is an auxiliary for fisquare_out
c
      t=tau0
      u2=u*u
      u3=u2*u
      t2=t*t
      t3=t2*t
      if (kl.eq.1) then
         dexu=dexp(-t/u)
         q=2d0*u3-u*dexu*(t2+2d0*t*u+2d0*u2)
      else if (kl.eq.2) then
         si2=si*si
         si3=si2*si
         if (si.ne.0d0) then
            dexs=dexp(-si*t)
            q=(2d0-dexs*(t2*si2+2d0*t*si+2d0))/si3
         else
            q=t3/3d0
         endif
      else if (kl.eq.3) then
         si2=si*si
         si3=si2*si
         if (si.ne.0d0) then
            dexs=dexp(-si*t)
            q=(si2*t2-2d0*si*t+2d0-2d0*dexs)/si3
         else
            q=t3/3d0
         endif
      else if (kl.eq.4) then
         dexu=dexp(-t/u)
         q=t2*u-2d0*t*u2+2d0*u3-2d0*u3*dexu
      else if (kl.eq.5) then
         dexu=dexp(-t/u)
         qx5=dexu*(6d0*u3+6d0*t*u2+3d0*t2*u+t3)-6d0*u3
         qt5=2d0*u2*dexu*(3d0*u+t)-6d0*u3+4d0*t*u2-t2*u
      endif
      return
      end



      subroutine qfun(q,tau0,r,t,s,a,b,nq)
      implicit real*8(a-h,o-z)
c     qfun finds the q-function at optical depths r and t
c     nq is the order of the Gaussian quadrature
c
      dimension s(nq),a(nq),b(nq)
      data eps/.95d-13/
      ii=1
      if (s(1).lt.eps) ii=2
      s1=0d0
      s2=0d0
      do 100 i=ii,nq
         be=s(i)
         et=dexp(-be*t)
         er=dexp(-be*r)
         ett=dexp(-be*(tau0-t))
         etr=dexp(-be*(tau0-r))
         be2=be*be
         s1=s1+a(i)*(er*(1d0+r*be)-et*(1d0+t*be))/be2
         s2=s2+b(i)*(ett*(t*be-1d0)-etr*(r*be-1d0))/be2
  100    continue
      if (s(1).gt.eps) then
         q1=0d0
         q2=0d0
      else
         r2=r*r
         t2=t*t
         r3=r2*r
         t3=t2*t
         q1=.5d0*a(1)*(t2-r2)
         q2=b(1)*(t3-r3)/3d0
      endif
      q=q1+q2+s1+s2
      return
      end
c
c
c
c
c
c
c
c        END OD HOMOGEN
c
c