30 August 1997 PFR, v 1.2, (c) G.J. Parker, 1997 USER MANUAL this document gives a rather brief overview on how to run the resonant radiation trapping program (PFR). documentation on every variable which can be set can be found in the accompanying reference manual. program control is determined by the variables which are set in the file "rad.in" which must reside in the directory from which the program is executed. unless otherwise noted, variables referred to below are in this data file. 1) geometry: a) slab_p, cyl_p, sph_p the simulation supports three geometries: infinite plane parallel, infinite cylindrical and spherical. the corresponding symmetries are assumed so that only one space coordinate (z, rho, r, respectively) is needed. the other two coordinates are automatically integrated over. infinite plane parallel: slab_p = .true., cyl_p = .false., sph_p = .false. infinite cylindrical: slab_p = .false., cyl_p = .true., sph_p = .false. spherical: slab_p = .false., cyl_p = .false., sph_p = .true. b) zmaxdp (cm) next the size of the system must be set. regardless of the geometry chosen, the size is set via the variable zmaxdp which is in cm. zmaxdp is the gap size for plane parallel or the outer radius for cylindrical/spherical geometry. 2) physical parameters: a) tvac (s) the vacuum lifetime of the resonant atoms is the time step of the simulation and therefore mostly irrelevant for the majority of the simulation since time is measured in units of the vacuum lifetime. b) k0 (cm^-1), k_lc_p the line center absorption coefficient is given by k0. k0 has been defined different ways by different authors. if k_lc_p is true the k0 is k(0), that is the value of the absorption coefficient at line center. since the lineshapes used in this simulations are normalized so that the integral over reduce frequency is 1, L(0) <> 1 for any lineshape L(x). k(x) = k'L(x) means that if the k0 is suppose to be k(0) (i.e. k_lc_p is true) then k0 = k' L(0) ==> k' = k0 / L(0). if k_lc_p is false, then k0 is k'. c) av for voigt or lorentz lineshape, av is the 'voigt/lorentz' parameter. it's not used for doppler profiles. in other words, for a lorentz profile, L_l(x) = (av/pi) / (x^2 + av^2) and for voigt, L_v(x) =(av / sqrt(pi^3)) * integral from -infty to +infty dw exp(-w^2) / ( (x-w)^2 + av^2). obviously, a doppler profile, L_d(x) = exp(-x^2)/sqrt(pi), is independent of av d) pcoll if this simulation is not a complete frequency redistribution (CFR) simulation, then the atom may or may not have a 'dephasing' collisions before it has a chance to re-emit the photon. pcoll gives the probability that such a collision will occur. 3) frequency redistribution. out of all the logicals in this section, only one should be set equal to true. all others should be set to false (w/ the possible exception of jw_dopp_p). for details on what happens if this is not true, see the REFERENCE MANUAL a) cfr_dopp_p, cfr_lrtz_p, cfr_vgt_p for CFR, one of the above logicals should be set corresponding to having CFR for a doppler (cfr_dopp_p), lorentz (cfr_lrtz_p) or voigt (cfr_vgt_p) lineshape, i.e. R(x, x') = L(x) L(x') b) jw_p (jw_dopp_p), pure_dopp_p, exact_coh_p, exact_incoh_p for PFR, current options are Jefferies-White approximation (jw_p AND not jw_dopp_p): 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; Jefferies-White but w/ a shifted doppler instead of delta (jw_p AND jw_dopp_p): R(x, x') = L_v(x') { (1-pcoll) { [1-a(x')] L_d(x) + a(x') L_d(x'-x)} + Pcoll L_v(x) }; pure doppler profile (pure_dopp_p): R(x, x') = 0.5 * erfc( x>), where erfc is the complementary error function and x> = max( x, x'); angle averaged frequency redistribution w/ coherent scattering (exact_coh_p): R(x,x') = (1-pcoll)/sqrt(pi^3) integral from (x> - x<)/2 to +infty du exp(-u^2) [atan( (x< + u)/av ) - atan( (x> - u)/av )] + pcoll * L_v(x), where x> = max( x, x'), x< = min( x, x'); angle averaged frequency redistribution w/ incoherent scattering (exact_incoh_p): R(x,x') = (1-pcoll)/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 )]+ pcoll * L_v(x). 4) type of simulation a) decay_p if the fundamental decay rate is wanted, then decay_p should .true. and prod_p .false. b) prod_p, prodfile if there is an external production and the steady-state resonant atom density profile is desired, then prod_p should be .true. while decay_p would be .false. prodfile is the name of a file (in the directory given by pathname if one is specified otherwise assumed to be in the current working directory) which gives the external production rate (1/cm^3/s) for each spatial cell. details of format of this data file can be found in the REFERENCE MANUAL. In a nutshell, the data file is scanned until at least four '>' at a start of a new line is encountered. then there must be izmax (see below) numbers. the first number is for the spatial cell closest to z (, rho or r) = 0 and the last number is for the spatial cell closest to z (rho or r) = zmaxdp. these numbers are the external production rate (1/cm^3/s) of resonant atoms. spatial cells are of uniform width. the rest of the data file is ignored. 5) computational mesh a) izmax the spatial mesh has cells of equal width and the number is determined by izmax. izmax must be less than the parameter zid (found in rad.inc) otherwise zid must be increased and the code recompiled. b) ikmax the frequency mesh has cells of UNEQUAL spacing which can either be set automatically or controlled by further user defined variables. ikmax is the total number of frequency cells and must be less than the parameter kid (found in rad.inc) otherwise kid must be increased and the code recompiled. c) freq_grid, ikbmax, freq_min, freq_mid, freq_max, k_min, k_mid, k_max the frequency mesh is constructed according to the value of freq_grid freq_grid = 0 generates the frequency mesh in an obsolete way- NOT recommended! freq_grid = 1 generates the 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) ) ). recommended for first (if not final) try. freq_grid = 3 generates the frequency mesh using the minimum frequency from line center (freq_min), frequency where high resolution should be around (freq_mid) and maximum frequency (freq_max). freq_grid = 2 generates the frequency mesh using the frequencies corresponding to the maximum k(x) (k_max), the k(x) where the finest resolution should be (k_mid) and the minimum k(x) (k_min). if ikbmax is < 2 or > ikmax-2, then ikbmax is ignored. otherwise ikbmax gives the number of frequency cells between line center and where the finest resolution should be. recommended to be set equal to -1. 6) internal simulation parameters a) ndt number of time steps (vacuum lifetimes) that will be performed to each call to the 'move' routine. a typical number would be 100, though 1 could be used just as well! output to the screen (for decay simulations- effective lifetime and number of resonant atoms left in discharge; for steady-state simulations- the number of resonant atoms in the discharge) is printed after each call to the mover. so, you want a decent number, otherwise you'll get a lot of output to the screen. b) print_ndt, filedat, esc_clear_p. after print_ndt calls to the mover (print_ndt * ndt time steps!), an output data file will be dumped in the directory given by pathname (current working directory if pathname is not defined) with name given filedat. each output will be appended with an ascending number. the output files give both the spatial profile of resonant atoms and the frequency spectrum of the escaping radiation. if esc_clear_p is true, then the escaping histogram spectrum is cleared after each printout. otherwise the spectrum is the sum through out the entire history of the run for each printout. c) num_ndt after num_ndt calls to the mover (num_ndt * ndt time steps), the code will stop. if this a decay simulation, the code would also stop if the number of resonant atoms drop to a level too small. d) n_phi_gauss, n_r_gauss if in cylindrical geometry, the spatial propagator must be found via numerical integration. Gaussian quadrature is chosen to do this. n_phi_gauss is the number of azimuthal points used in the integration while n_r_gauss is the number used in the radial integration. the maximum values are 300 and 20, respectively. recommended values would be 300 and 1. e) nave in the calculation of the spatial propagator, nave gives the number of points that the INITIAL cell should be averaged over. for slab and spherical geometries, this can be large (~20 works well) and it's not expensive. for cylindrical, it should be at least 2 but large values will take quite awhile to compute the propagator. if you have time or a fast computer, large values are typically better. (don't know since i have neither!) f) read_p, oldfile instead of recalculating the propagators, you can read in a previously saved one. this is done if read_p is true and oldfile contains a name of such a data file (in the directory pathname or in the current working directory if pathname is not set). if the file does exist, then the following variables are reset to the values specified in the data file oldfile: zmaxdp, k0, tvac, izmax, nave, slab_p, cyl_p, sph_p, ikmax, ikbmax, av, pcoll, jw_p, jw_dopp_p, pure_dopp_p, cfr_dopp_p, cfr_lrtz_p, cfr_vgt_p, exact_coh_p, exact_incoh_p, n_phi_gauss, n_r_gauss g) save_p, newfile if save_p is true, then the propagators will be saved to the data file newfile (in the directory given by pathname or in the current working directory if not set). 7) Analysis at this time, you should be able to execute the code and hopefully get some output. the question is then whether you believe what the code is saying. a) decay simulations: during the simulation, the code will print out the instantaneous decay rate. this number should go smoothly to some fixed value and stay there. if it doesn't, you should increase the duration of the simulation (num_ndt) and try again. if it still doesn't settle down, you may want to change the frequency mesh (see below) b) steady-state simulations: during the simulation, the code will print out the instaneous total number of resonant atoms in the discharge. this number should eventually reach a fixed number. if not, increase num_ndt until it does. once you are reasonably convinced that the code has found a steady-state solution or the fundamental decay rate, looking at the last output file (filedat) is essential. looking at the escaping frequency spectrum, the peak (the ratio column) should be close to the finest resolution in frequency and drop off smoothly on each side. at large frequency, this number should be close to zero. if it isn't, you may want to increase the maximum frequency (or similarly decrease the minimum k(x)) by using freq_grid. on the other hand, if there are many large frequency cells which are zero (or *very* small), then you may want to decrease (increase) the maximum frequency (minimum k(x)). near line center, if many cells are zero and this is not an extremely opaque case, you may want to increase (decrease) the minimum frequency (maximum k(x)). EOF