PHOTON USE
   AUTOPLOT
   file
   PHI 5
 
   cl;d 1 h1;d 1 ha;col3 1;blb4 2;redr
   msg    temperature profile; press  to continue
   pause
   cl;d 1 radx;d 1 ra;col3 1;blb4 2;redr
   msg    radiation-flux profile; press e to end
   pause;end
   ENDUSE
    GROUP 1. Run title and other preliminaries
TEXT(X-D Radiative Equilibrium In Slab 
TITLE
  DISPLAY
   The problem considered is that of radiative heat transfer
   in a 1d plane-parallel slab in radiative equilibrium. The
   slab is of thickness L and the gray medium may absorb,
   emit and scatter radiation. At x=L the wall is a diffuse
   emitter and reflecter kept at a fixed temperature. At the x=0
   the net radiative heat flux Qrad is specified at the wall.
  ENDDIS
 
   Since the energy transfer is by pure radiation, the energy
   equation is given by:
 
        -d/dx (Qrad) = 0
 
   The equation for the composite radiative heat flux is given
   by:
 
        d/dx ( 1/(a+s) d/dx (Rx) ) + a (E - Rx) = 0
 
   where a is the absorption coefficient, s is the scattering
   coefficient, Rx is the composite radiation flux defined as
   the average of the +ve and -ve radiation fluxes, and E is
   the black-body emissive power. It may be noted that the
   radiative heat flux is given by:
 
        d/dx (Qrad) = 2a (E - Rx)
 
   The black-body emissive power E=sig*T**4  where sig is the
   Stefan-Boltzmann constant and T is the temperature of the
   medium. The problem is the determination of the temperature
   and composite radiative-flux distributions, as given by the
   following analytical solutions:
 
     Rx = E   ( radiative equilibrium )
 
     E  = Ew + Qrad*[ 1./emw - 0.5 + 0.5*(a+s)*L*(1 - x/L) ]
 
    where Ew is the emissive power at the wall and emw is the
    emissivity of the wall. For the case considered below, Qrad
    is taken to be Qrad = 2.*Ew.
 
    The locally-defined parameters are as follows:
        GSIGMA     Stefan-Boltzmann constant       {W/m**2/K**4}
        SCAT       Scattering coefficient          { 1/m  }
        ABSORB     Absorption coefficient          { 1/m  }
        EMIW       emissivity of the wall
        TWAL       wall temperature at x=L         {  K   }
        QRAD       net radiative heat flux at x=0  {W/m**2}
 
CHAR(CH1);REAL(GSIGMA,SCAT,ABSORB,EMIW,TWAL,QRAD)
GSIGMA=5.6697E-8;SCAT=0.5;ABSORB=0.5;EMIW=1.0;TWAL=1000.0
QRAD=2.*GSIGMA*TWAL**4
    GROUP 3,4,5. X,Y,Z-direction grid specification
GRDPWR(X,50,1.0,1.0)
    GROUP 7. Variables stored, solved & named
CP1=1.0
MESG( Enter required energy variable ? (TEM1 or H1)
IF(:CH1:.EQ.) THEN
+  READVDU(CH1,CHAR,H1)
ENDIF  
IF(:CH1:.EQ.TEM1) THEN
+ MESG( TEM1 solution selected
ELSE
+ MESG( H1 solution selected
+ TMP1=LINH;TMP1B=1.0/CP1
ENDIF
RADIAT(FLUX,ABSORB,SCAT,:CH1:);STORE(EMPO)
    GROUP 8. Terms (in differential equations) & devices
   ** Deactive conduction & any built-in sources
TERMS(:CH1:,N,N,N,N,P,P)
    GROUP 11. Initialization of variable or porosity fields
FIINIT(RADX)=0.5*QRAD;FIINIT(:CH1:)=TWAL
  ** analytical solution
REAL(EW,EG,QRADA,ALF,BET,TA,RAN,GX);STORE(HA,RA);INTEGER(JJM1)
EW=GSIGMA*TWAL**4;QRADA=2.*EW;ALF=1./EMIW-0.5
BET=0.5*(ABSORB+SCAT)*XULAST
DO JJ=1,NX
+PATCH(IN:JJ:,INIVAL,JJ,JJ,1,NY,1,NZ,1,1)
+GX=0.5*XFRAC(JJ)
IF(JJ.NE.1) THEN
+JJM1=JJ-1;GX=XFRAC(JJM1)+0.5*(XFRAC(JJ)-XFRAC(JJM1))
ENDIF
+GX=GX*XULAST;EG=EW+QRADA*(ALF+BET*(1.-(GX/XULAST)))
+RAN=EG;TA=(EG/GSIGMA)**0.25
+INIT(IN:JJ:,RA,ZERO,RAN);INIT(IN:JJ:,HA,ZERO,TA)
ENDDO
    GROUP 13. Boundary conditions and special sources
   ** Net radiation flux from wall
PATCH(WALLRA,WEST,1,1,1,NY,1,NZ,1,1)
COVAL(WALLRA,RADX,FIXFLU,0.5*QRAD)
PATCH(WALLRB,EAST,NX,NX,1,NY,1,NZ,1,1)
COVAL(WALLRB,RADX,EMIW/(2.0-EMIW),GSIGMA*TWAL**4)
    GROUP 15. Termination of sweeps
LSWEEP=50
    GROUP 16. Termination of iterations
RESREF(:CH1:)=1.E-12*QRAD;RESREF(RADX)=0.5*RESREF(:CH1:)
    GROUP 17. Under-relaxation devices
RELAX(:CH1:,FALSDT,100./QRAD)
    GROUP 22. Spot-value print-out
IXMON=NX/2;NPLT=5;NXPRIN=NX/10
    GROUP 23. Field print-out and plot control
OUTPUT(:CH1:,Y,N,N,Y,Y,Y)
PATCH(XWISE,PROFIL,1,NX,1,1,1,1,1,1)
PLOT(XWISE,:CH1:,0.0,0.0);PLOT(XWISE,RADX,0.0,0.0)
    GROUP 24. Dumps for restarts
TSTSWP=-1