#immersol)
#declare
   IMMERSOL Emitting and absorbing gray media
            (1D, TEM1)
TEXT(IMMERSOL 1D Radiative equilibrium   :202
TITLE
  **************************************************************
   PHOTON USE
   AUTOPLOT
   FILE
   PHI 5
 
   CL;DA 1 TEM1;DA 1 T3;COL3 1;BLB4 2;REDR
   MSG    Temperature (curve) and T3 (crosses) profiles
   ENDUSE
  **************************************************************
#cls
  DISPLAY
   The problem is as for case 201; but now allowance is made for
   absorption and emission by the medium between the two plates.
   The gap between the plates is WG (default = 1.0), and the
   absorption/emission coefficient of the medium is EMISS.
 
   The exact solution for the heat flux is:
         QRAD = (EH - EC)/(0.75*ROSS*WGAP + 1/EMISH + 1/EMISC - 1)
   where EH and EC are the black-body emissive powers of the hot and 
         cold surfaces, i.e. temperature**4 * Stefan-Boltzman 
         constant;
         ROSS=(EMISS+SCATT) is the Rosseland coefficient.
         The ROSS*WG product gives optical thickness of the gas slab
 
  ENDDIS
   The case is set to exclude influence of the heat conduction in
   a gas. The user can activate heat conduction by substituting
   command TERMS(TEM1,N,N,N,P,P,P) by TERMS(TEM1,N,N,Y,P,P,P) .
#direct
#temps
MESG(  ENTER optical thickness Kr*Wg (default 1.):
READVDU(OPTHI,REAL,1.0)
  settings of wall thicknesses and the gap between them
WSL1 = 0.01;  WSL2= 0.01;  WGAP= 1.0;  
  conductivity of medium swet to zero
KGAS= 0.0
  QRAD first used to store EH-EC (see above)
QRAD= SIGMA*(THOT**4 - TCLD**4)
  KROSS deduced from above-set optical thickness; scattering
  coefficient set to zero; then medium emissivity deduced
KROSS= OPTHI/WGAP;  SCATT= 0.0;  EMISS= KROSS-SCATT
MESG( Expected radiative heat flux (W/m**2):
  QRAD now computed, based on above settings.
  Note that KGAS has been set to zero, so that NN1 and NN2 = 0
NN1 = KGAS*KROSS/4./SIGMA/THOT**3
NN2 = KGAS*KROSS/4./SIGMA/TCLD**3
QRAD= QRAD + KGAS/WGAP*(THOT-TCLD)*0.75*OPTHI
AA1 = (1./EMISH-0.5)/(1.+0.75*NN1)
AA2 = (1./EMISC-0.5)/(1.+0.75*NN2)
QRAD= QRAD/(0.75*OPTHI + AA1 + AA2)
QRAD
MESG( Expected gas temperatures near walls are:
TGCL = ((SIGMA*TCLD**4 + QRAD*AA2)/SIGMA)**0.25
TGHT = ((SIGMA*THOT**4 - QRAD*AA1)/SIGMA)**0.25
TGCL
TGHT
#pause
SOLVE(TEM1)
STORE(PRPS)
#geom1
#radflux
TERMS(TEM1,N,N,N,N,Y,N);  PRNDTL(TEM1)=1.e10
INIADD= F;  FIINIT(TEM1)= (THOT+TCLD)/2.;  FIINIT(PRPS)= 0.
FIINIT(T3)=FIINIT(TEM1)
 COVAL(HOT,TEM1,FIXVAL,THOT); COVAL(COLD,TEM1,FIXVAL,TCLD)
 INIT(SOL1,PRPS,0.0,111.); INIT(SOL2,PRPS,0.0,112.)
LSWEEP= 500
IXMON=NX/2+1;IYMON=NY/2+1;IZMON=NZ/2+1
NXPRIN=1; NYPRIN=1; NZPRIN=1
OUTPUT(LTLS,N,N,N,N,N,N); OUTPUT(WDIS,N,N,N,N,N,N)
  
  load macro for storing and setting emissivity and scattering 
  coefficient, using EMISS, SCATT, EMISH and EMISC
load($r193)

VARMAX(T3)=THOT;VARMIN(T3)=TCLD
VARMAX(TEM1)=THOT;VARMIN(TEM1)=TCLD
   Activating the following In-Form statements 
   (by moving them 2 spaces to the left) 
   enables the coefficients in the equations for TEM1 and T3 
   to be inspected.
   
    inform7begin
  (stored var #3-1 is t3-tem1)
  (stored var #rs1 is resi(tem1))
  (stored var #rs3 is resi(t3))
  (stored var #ap1 is apco(tem1))
  (stored var #ap3 is apco(t3))
  if(idir.eq.1) then
  (stored var #ae1 is aeco(tem1))
  (stored var #ae3 is aeco(t3))
  endif 
  if(idir.eq.2) then
  (stored var #an1 is anco(tem1))
  (stored var #an3 is anco(t3))
  endif 
  if(idir.eq.3) then
  (stored var #ah1 is ahco(tem1))
  (stored var #ah3 is ahco(t3))
  endif 
    inform7end