23 August 1997 18 July 1997 04 July 1997 02 July 1997 10 January 1997 PFR, v 1.2, (c) G.J. Parker, 1997 REFERENCE MANUAL this file contains documentation for the propagator method for resonant radiation trapping with partial (and complete) frequency redistribution. the program will be referred to as PFR the accompaning file title USER MANUAL describes how to run this program while this file documents all of the variables which can be set in the input file. this program is based, in part, on the following papers: J.E. Lawler, G.J. Parker and W.N.G. Hitchon, "Radiation trapping simulations using the propagator function method", J. Quant. Spectros. Radiat. Transfer. v49, p627 (1993). G.J. Parker, W.N.G. Hitchon and J.E. Lawler, "Radiation trapping simulations using the propagator function method: complete and partial frequency redistribution", J. Phys. B v26, p4643, (1993). A.F. Molisch, G.J. Parker, B.P. Oehry, W. Schupita and G. Magerl, "Radiation trapping with partial frequency redistribution: comparison of approximations and exact solutions", J. Quant. Spectros. Radiat. Transfer v53, p269, (1995). the PFR is free, however it must not be given to anyone with out my knowledge and consent. currently, i can be reached at parker9@llnl.gov and/or parker@percy.engr.wisc.edu. acknowledgements in any work that uses this program in whole or in part is also required. other than these caveats, you can do what you wish with this code. source files: ------------ rad.f- driver routine (fortran source) inputp.f- set up routine for PFR (fortran source) prob.f- routine to calculate propagators J and Q (fortran source) mov.f- routine to step density by a lifetime (fortran source) misc.f - misc. routines for output and functions (fortran source) rad.inc- include file for common blocks (fortran source) Makefile- make file to compile program- works for HP (make file) input files: ----------- rad.in- namelist file to define mesh and other parameters and numerous other files which are specified in the above file. revision history: ---------------- 0.1 1993 - 1996 numerous partial releases 1.0 10 Jan. 1997 first release 1.0.1 2 July 1997 added ability to adjust number of Gaussian quadrature points. 1.0.2 4 July 1997 enhanced output files and screen output. made the determination of frequency mesh automatic on user description. 1.1 18 July 1997 enhanced user specification of frequency mesh 1.2 23 Aug. 1997 fixed frequency redistribution propagator for exact (coherent and incoherent scattering). Also possible fix for JW plus shifted Doppler for extremely transparent gases philosophy of code: ------------------ the PFR is based on a series of papers published while i was at the University of Wisconsin- Madison and later while i was at LLNL. the original paper used a different code in that complete frequency redistribution was assumed with a Lortenzian lineshape in slab (plane parallel) and spherical geometries. quickly after that, my colleagues and i moved onto partial frequency redistribution using the standard Jefferies-White approximation and a Voigt lineshape and we included cylindrical geometry. after moving onto LLNL, on the encouragement of A.F. Molisch, i revived the code and incorporated the angle averaged exact redistribution for a pure Doppler lineshape and for a Voigt lineshape. since then, i have extended the program to also do CFR also. the program is generally efficient but because of the lack of analytic results for cylindrical geometry, i was forced to use numerical integration (ie Gaussian quadrature) to do the spatial integration to find the Q propagator, accordingly it can take quite awhile to compute. additionally, the exact angle averaged coherent scattering frequency redistribution propagator also must be integrated numerically and can be expensive. other than these two cases, construction of the propagators is generally cheap and the majority of the cpu time is spent finding either the steady-state density profile (if there is an external source) or the decay rate. this is simply a matrix multiplication and could be optimized for the specific computer at hand (see mov.f). definition of problem: --------------------- The PFR solves the generalized Holstein-Biberman equation for resonant radiation trapping in slab, cylindrical or spherical geometries with Doppler, Voigt, or Lorentz lineshape in the Jefferies-White, modified Jefferies-White, exact angle averaged (either coherent or incoherent scattering) (excluding Lorentzian lineshape) or CFR approximations. an arbitrary external source of resonant atoms can be specified and the steady state density profile can be solved for or the ratio of trapped decay rate to that of vacuum decay rate can be found. the lineshape is assumed to by symmetric and a simple standard shape- either Doppler, Lorentzian or Voigt. the method of solution is to cast the Holstein-Biberman equation in terms of kernels or propagators of the intregal-differential equation and then vectorize phase space (space and frequency) which leads to a set of coupled algebraic equations for the resonant atom densities and frequency spectrum. generation and execution of executable: -------------------------------------- typing either 'make rad' or 'make rad_d' at the command line prompt in the directory which contains the six source files and the Makefile will generate executable code. the latter command will produce a slower but debuggable version of the code. to execute the code, the current directory must contain the input file rad.in. the executable code need not be in the same directory. input file definitions: ----------------------- as mentioned previously, there is one main input file for the PFR: rad.in this file defines the geometry, parameters of the problem, the manner in which the PFR is executed. other input/output files (e.g. external production file) are defined and subsequently read. >>>> rad.in <<<< this file is the first data file read and must exist. the valid entries in this file are: ****** zmaxdp = radius/gap size (cm) tvac = vacuum lifetime of excited resonant atom (s) k_lc_p = .true. --> k0 below is the value of k at x = 0 .false. --> k0 below is the absorption coefficient integrated over reduced frequency, ie k(x) = k0 L(x) k0 = absorption coefficient (cm^-1) av = Voigt parameter (#) Pcoll = probability of a dephasing collision before emission (#) jw_p = use Jefferies-White approximation? jw_dopp_p = if jw_p is true, modify JW by replacing delta function w/ shifted Doppler? pure_dopp_p = use Doppler profile for line shape- NOT Voigt ? cfr_dopp_p = CFR with a Doppler profile? ie L_d(x) = exp(-x^2)/sqrt(pi) cfr_lrtz_p = CFR with a Lorentzian profile? ie L_l(x) = (av/pi) / (x^2 + av^2) cfr_vgt_p = CFR with a Voigt profile? ie L_v(x) = (av / sqrt(pi^3)) * integral from -infty to +infty dw exp(-w^2) / ( (x-w)^2 + av^2) exact_coh_p = use angle averaged frequency redistribution w/ coherent scattering? exact_incoh_p = use angle averaged frequency redistribution w/ incoherent scattering? izmax = number of spatial cells (#) ikmax = number of frequency cells (#) ikbmax = number of frequency cells between line center and x* where k(x*) zmaxdp = 1 or user specified(#) if less than 1, ikbmax is determined automatically by code freq_grid = integer from 0 to 3- specifies how frequency mesh should be generated (#) freq_min = minimum reduced frequency freq_mid = mid-range reduced frequency freq_max = maximum reduced frequency k_min = minimum value of k k_mid = mid-range value of k k_max = maximum value of k nave = initial cell is broken into this many subcells when computing the spatial propagator Q (#) num_ndt = this number times ndt is the maximum # of time steps (#) print_ndt = this number times ndt is the number of time steps between printouts (#) ndt = number of time steps for each call to the 'mover' (#) esc_clear_p = zero out escaping spectra histogram after each printout? decay_p = is this a 'decay' simulation- ie no external source? prod_p = is this a 'steady-state' simulation? is there an external source? slab_p = is this a slab geometry? cyl_p = is this a cylindrical geometry? sph_p = is this a spherical geometry? n_phi_gauss = number of azimuthal points in Gaussian quadrature n_r_gauss = number of radial points in Gaussian quadrature pathname = path to where input/output files are prodfile = name of production file if this is a prod_p run filedat = output file root name read_p = read in a previously saved J and Q propagators? save_p = save J and Q propagators? oldfile = if read_p is true, file where J and Q are stored newfile = if save_p is true, file where J and Q should be saved ****** zmaxdp is the radius or gap size in cm. default value is 1.0 cm tvac is the vacuum lifetime of excited resonant atom (s). default value is 1.0 s k_lc_p tells how k0 (see below) is defined. if true then k0 is the value of the absorption coefficient at line center, ie k0 = k(0) = k' L(0), where L(x) is specified lineshape. if false, then k0 is k' in this expression. default is false. k0 is the absorption coefficient either at line center or integrated over reduced frequency. which one is determined by k_lc_p (see above) in cm^-1. default is 1.e5 cm^-1 av is the Voigt parameter (#). default value is 0.1 Pcoll is probability of a dephasing collision before emission. this is ignore if the run is a CFR run (ie if any of the following is true: cfr_dopp_p, cfr_lrtz_p or cfr_vgt_p) or if a pure Doppler lineshape is used (ie pure_dopp_p is true). default value is 0.5 the following logicals are mutually exclusive (ie only one of them is true): jw_p, pure_dopp_p, cfr_dopp_p, cfr_lrtz_p, cfr_vgt_p, exact_coh_p and exact_coh_p. if there is more than one which is true or none of them is true the following logic, in the order given, is followed: if jw_p is true, then all others are set to false if exact_coh_p is true, then all others are set to false if exact_incoh_p is true, then all others are set to false if pure_dopp_p is true, then all others are set to false if cfr_dopp_p is true, then all others are set to false if cfr_lrtz_p is true, then all others are set to false if cfr_vgt_p is true, then all others are set to false otherwise all set to false except for jw_p which is set to true jw_p is true if the Jefferies-White approximation should be used. a Voigt lineshape is assumed. R(x, x') = L_v(x') { (1-pcoll) { [1-a(x')] L_d(x) + a(x') delta(x'-x)} + Pcoll L_v(x) }, where a(x) = max(0, min(1, 1 - L_d(x) / L_v(x) )) and delta(x) is dirac delta function. default is false. jw_dopp_p is true if Jefferies-White approximation should be modified by replacing the delta function w/ shifted Doppler lineshape. only used if jw_p is true. R(x, x') = L_v(x') { (1-pcoll) { [1-a(x')] L_d(x) + a(x') L_d(x'-x)} + Pcoll L_v(x) } default is false. pure_dopp_p is true if a pure Doppler profile should be used. Pcoll, av is ignored. R(x, x') = 0.5 * erfc( x>), where erfc is the complementary error function and x> = max( x, x'). default is false. cfr_dopp_p is true if CFR with a Doppler profile should be used. Pcoll and av ignored. R(x', x') = L_d(x) L_d(x'). default is false. cfr_lrtz_p is true if CFR with a Lorentzian profile should be used. pcoll and av ignored. R(x', x') = L_l(x) L_l(x'). default is false. cfr_vgt_p is true if CFR with a Voigt profile should be used. pcoll ignored. R(x', x') = L_v(x) L_v(x'). default is false. exact_coh_p is true if angle averaged frequency redistribution w/ coherent scattering should be used. R(x,x') = 1/sqrt(pi^3) integral from (x> - x<)/2 to +infty du exp(-u^2) [atan( (x< + u)/av ) - atan( (x> - u)/av )], where x> = max( x, x'), x< = min( x, x'). default is false. exact_incoh_p is true if angle averaged frequency redistribution w/ incoherent scattering should be used. R(x,x') = 1/sqrt(pi^5) integral from 0 to +infty du exp(-u^2) [atan( (x + u)/av ) - atan( (x - u)/av )] [atan( (x' + u)/av ) - atan( (x' - u)/av )]. default is false. izmax is the number of spatial cells (#). default is the maximum array is declared: zid, see rad.inc. ikmax is the number of frequency cells (#). default is the maximum array is declared: kid, see rad.inc. ikbmax is the number of frequency cells between line center and x* where k(x*) zmaxdp = 1 (#) or freq_mid or k_mid = k(x*), see below. if ikbmax < 1, then ikbmax is determined automatically. default is -1. freq_grid is an integer between 0 and 3, inclusive. if outside of this range, it's set to 1. the values mean the following: 0 -> generate frequency mesh old way using ikbmax- NOT recommended 1 -> generate frequency mesh automatically based on variation of k(x). this is the same as freq_grid = 2 with k_max = min( k0, -2 ln(1.e-6) / (zmaxdp/izmax) ), k_min = - k0 ln(0.99) / (2 zmaxdp), k_mid = exp( 0.5 ( ln(k_min) + ln(k_max) ) ) 3 -> use the frequencies specified by freq_min, freq_mid and freq_max for the minimum frequency, frequency where high resolution should be around and maximum frequency. if any of these values less than 0, then use value computed as in freq_grid=1. 2 -> use the frequencies corresponding to k(x) = k_max, k_mid, k_min. if any of these values less than 0, then use value computed as in freq_grid=1. the frequency cells will be the finest around freq_mid/k_mid. default is 1 freq_min is minimum reduced frequency if freq_grid = 3. default is -1. freq_mid is mid-range reduced frequency if freq_grid = 3. default is -1. freq_max is maximum reduced frequency if freq_grid = 3. default is -1. k_min is minimum value of k(x) if freq_grid = 2. default is -1. k_mid is mid-range value of k(x) if freq_grid = 2. default is -1. k_max is maximum value of k(x) if freq_grid = 2. default is -1. nave is the number of subcells each initial cell is broken into when computing the spatial propagator Q (#). default is 2. num_ndt is the maximum calls to the 'mover'. default is 10000. print_ndt is the number times mover is called between printouts. default is 10. ndt is number of time steps (vacuum decay lifetimes) for each call to the 'mover'. default is 100. esc_clear_p is true if escaping spectra histogram should be reset after each printout. default is true. decay_p and prod_p are mutually exclusive. prod_p is true if steady-state calculation is to be performed. if so the prodfile must exist in the directory given by pathname. if it doesn't, prod_p is set to false and decay_p is set to true. otherwise, decay_p is set to false. default is false. decay_p is true if this a 'decay' simulation. default is true. slab_p, cyl_p and sph_p are mutually exclusive. if all are false or more than one is true the following logic is used: if cyl_p is true, then others are set to false. if slab_p is true, then others are set to false. if sph_p is true, then others are set to false. otherwise cyl_p is set to true and others to false. slab_p is true for slab geometry. default is false. cyl_p is true for cylindrical geometry. default is true. sph_p is true for spherical geometry. default is false. n_phi_gauss is the number of azimuthal points in Gaussian quadrature if the geometry is cylindrical (ie cyl_p = true), otherwise ignored. default is 300, range between 50 and 400. n_r_gauss is number of radial points in Gaussian quadrature if the geometry is cylindrical (ie cyl_p = true), otherwise ignored. default is max(1, int(nave / 2) ). pathname is the directory for prodfile and output files. if set to null then current directory is used. default is null. prodfile is the name of production file if prod_p is true. it's in the directory given by pathname. default is 'none'. filedat is the name of output files. default is 'rad.out'. read_p is true if a previously stored J and Q propagators are to be read in. if true, file is in the directory given by pathname and file is given in oldfile. if there is an error trying to read this file, read_p is set to false and the run continues. otherwise all of the variables in this file are IGNORED (ie set to the values which are stored in this data file) EXCEPT for the following list: num_ndt, print_ndt, ndt, decay_p, prod_p, prodfile, filedat, pathname, esc_clear_p, read_p, save_p, oldfile and newfile. default is false save_p is true if J and Q propagators are to be saved. if true, data will be saved to the file specified by newfile in the directory given by pathname. default is true. oldfile is the data file in the directory given by pathname to be read if read_p is true. default is 'oldprop.dat'. newfile is the data file in the directory given by pathname where J and Q are to be saved if save_p is true. default is 'newprop.dat'. >>>> prodfile <<<< this string is the name of the external production rate of resonant atoms. this file is required if prod_p is true. if it doesn't exist in the directory given by pathname, prod_p is set to false and decay_p is set to true. if the file exists, it is scanned for valid entries. a valid entry starts with at least four (4) '>' on a newline. on the next non-blank line should be the production rate (1/cm^3/s) of the external source for the leftmost cell in slab geometry or the inner most radius cell for cylindrical/spherical geometry. there should be izmax-1 more numbers for each remaining cell moving to larger z or r, respectively. there must be at least izmax numbers. anything after the last number is ignored. output: ------ PFR writes to standard output device what it has read from rad.in and/or oldfile as the case may be. it then proceeds to the calculation of J and Q if they are not read in. if prod_p is true and prodfile was successfully read in, the program will then produce output files according to ndt, print_ndt and num_ndt and will exit once the specified number of time steps have been executed. if decay_p is true, then after each call to the 'mover', the current ratio of decay rate to vacuum decay rate is printed along with the total number of excited resonant atoms left in the volume. also output files will be saved according to ndt, print_ndt and num_ndt. the simulation will continue until the number of specified time steps have been executed or the number of excited resonant atoms in the discharge falls below 1.e-22. [this number is not an integer!] output files are named according to filedat with an integer appended to the file name for each printout in ascending order starting at 0. the first section of the output files show the density profile of excited resonant atoms a a function of position. density, number and density / maximum density are all shown as a function of position (z or r). densities are in 1/cm^3 the section section shows the spectrum escaping the system. the density is simply the number of photons divided by the frequency cell width. number is the raw number of photons in any given frequency cell. number is not shown for slab geometry, but spectrum density for each plane IS shown. if esc_clear_p is true, then this spectrum is cleared after each printout. otherwise it's the running sum of escaping photons. bug reports: ----------- known bugs: 1) for small values of k0, the spatial propagator Q may become numerically unstable. this manifests itself by giving incorrect density profiles and decay rates. this is simply a problem of the accuracy in the expressions needed in computing Q. more accurate expressions could result in better results. 2) excessive values of nave can lead to strange results. safe value is 2. when you run across a bug/error (notice i said when), please fill out a brief but informative report and send to me. also, if you have squashed a bug which is in the original code, i would like to know that also. EOF