----------------------------------
+
+   deriving equations 32-34.
+
----------------------------------

1. I am able to rederive eqn (32) if I make the following assumptions:

   a. The upward monochomatic intensity at the lower bountry,
      I_nu(mu, z0), is given by I_nu(mu,z0) = B_nu(z0).

   b. The downward intensity is ignored (hence the term, "cool to space").

   c. The left-hand side of eqn 33 is B(z0), not B.


   Start off with the monochromatic equation for the source function,
assuming no cascading.

  S_nu = B_nu * (J/B + e)/(1+e)                                 (RVY-7)
  
Multiply both sides by (k_nu/k0) and integrate over frequency, and use 
eqs (RVY-6) and (RVY-7) to get

  S = B * (J/B + e)/(1+e) = (J + e B)/(1+e)                         (1) 

Now we turn to equation 32.  The second term of eqn 1, e B/(1+e), 
is the first term of eqn (RVY-32).  Now we just need an expression for J.

  Formal integration of the plane-parallel RT equation (RVY-31)
gives the monochromatic intensity as a function of mu, the cosine of
the zenith angle.


  for 0 < mu < 1

  I_nu(tau, mu) = I(tau0) exp(-(tau0-tau)/mu) 

                   / tau0
                 + |       S_nu(tau', mu) exp(-(tau'-tau)/mu) dtau'/mu
                   / tau

  for -1 < mu < 0
                   / 0
  I_nu(tau, mu) =  |       S_nu(tau', mu) exp(-(tau'-tau)/mu) dtau'/mu
                   / tau


Make the following substitutions:

   S_mu(tau', mu) is independent of mu.

   I(tau0) = B(z0).  The lower boundary is a blackbody.

   dtau' = - k_nu(z') dz'

   tau(z,z') = | tau' - tau |

and get


  for 0 < mu < 1

  I_nu(tau, mu) = B(z0) exp(-tau(z0,z)/mu) 

                   / z
                 + |       S_nu(z') exp(-tau(z,z')/mu) k_nu dz'/mu
                   / z0

  for -1 < mu < 0
                    / inf
  I_nu(tau, mu) = - |       S_nu(z') exp(+tau(z',z)/mu) dtau'/mu
                    / z


So, we see that 

            1  /  0                    1  / +1  
  J_nu(z) = -  |     I_nu(z, mu) dmu + -  |     I_nu(z, mu) dmu
            2  / -1                    2  /  0

has 3 terms.  Let's call them J1_nu, J2_nu, and J3_nu.


                     / +1
    J1_nu(z) = (1/2) |     B(z0) exp(-tau(z0,z)/mu) dmu
                     / 0


                     / +1      / z
    J2_nu(z) = (1/2) |    dmu  |       S_nu(z') exp(-tau(z,z')/mu) k_nu dz'/mu
                     / 0       / z0


                     /  0      / inf
   J3_nu(z) = -(1/2) |    dmu  |       S_nu(z') exp(+tau(z',z)/mu) k_nu dz'/mu
                     / -1      / z


Take each one in turn, remembering the definitions for the exponential
integrals from Abramowitz and Stegun (5.1.4)

              / inf
    En(tau) = |     u^-n exp(-tau u) du
              / 1


Using the substitution 1/mu -> u, we derive

  
              / 1
    E1(tau) = |     (1/mu) exp(-tau/mu) dmu
              / 0

  
              / 1
    E2(tau) = |     exp(-tau/mu) dmu
              / 0
 

Now we evaluate J1_nu, pulling B(z0) out of the integral.
                     
    J1_nu(z) = (1/2) B(z0) E2( tau(z0,z) )
                     
This is just eqn (RVY-33), which must also be B0 in eqn (RVY-32).


Next we evaluate J2_nu, rearranging the order of the integral

                     / z                     / +1
    J2_nu(z) = (1/2) |    dz 'S_nu(z') k_nu  |    dmu exp(-tau(z,z')/mu)/mu
                     / z0                    / 0

and we evaluate the inner integral

                     / z                     
    J2_nu(z) = (1/2) |    dz' S_nu(z') k_nu  E1( tau(z,z') )
                     / z0                    

To evaluate J3_nu, we also switch the order of the integrals, and also
make the substitution mu -> -mu

                     /  inf                     / 1
    J3_nu(z) = (1/2) |      dz' S_nu(z')  k_nu  |   dmu  exp(-tau(z',z)/mu)/mu
                     / z                        / 0

giving, again,

                     / inf                     
    J3_nu(z) = (1/2) |    dz' S_nu(z') k_nu  E1( tau(z,z') )
                     / z                    


J2_nu and J3_nu can be combined, bearing in mind that tau(z,z1) will
be calculated differently for z > z' and z < z'.  Let's call it J4_nu

                     / inf                     
    J4_nu(z) = (1/2) |    dz' S_nu(z') k_nu  E1( tau(z,z') )
                     / z0                    


To turn the monochromatic mean intensity into the mean intensity of the
band, we multiply by k_nu(z), integrate over frequency, and divide by k0(z).