To read the mathematics in this version,
which was produced by
TtH,
under Unix and Netscape,
put the following line in .Xdefaults or .Xresources:
Netscape*documentFonts.charset*adobefontspecific: iso88591
Nov 19, 2000
Lower hybrid heating has been of interest for some time. Computational analysis perhaps started with an approach to calculating the spectrum of waves launched from a coupler [1]. The eikonal method [2], or geometrical optics [3], also referred to as ray tracing, was known to be available to problems in plasma physics, and several workers started investigating the consequenses of the ray method [4,5,6,7,8]. Success at driving the toroidal current in a tokamak [9,10] following theoretical predictions [11] caused increased interest in theory [12,13], modelling [14,15,16,17] and related considerations [18,19,20,21,22]. A natural extension was to add to the comprehensive analysis tools such as TRANSP [23,24] and TSC [25] the the best models of lower hybrid heating (LHH) and lower hybrid current drive (LHCD). Such an extension was particularly desirable becasue of a possible investment in a major new project to operate a tokamak in a steady state mode [26,27].
A static but magnetically accurate model [28,29] plus studies of wave chaos and accessibility [18,22,30,31] revealed many of the important issues.
The Lower hybrid Simulation Code (LSC) is a computational model of lower hybrid current drive in the presence of an electric field [32] derived from direct forerunners in references [8] and [15]. Heating of electrons and ions is also computed. Details of geometry and the plasma profiles are treated, so that the accessibility problems are treated as fully as possible in the ray model [32,21]. Twodimensional velocity space effects are approximated in a onedimensional FokkerPlanck treatment [12]. LSC contains a model of an xray camera to help compare results of the model with the physical camera [33,34,35,36].
Applications include work in references [32], [37] and [38].
The LSC was originally written to be a module for lower hybrid current drive and heating called by the Tokamak Simulation Code (TSC) , which is a numerical model of an axisymmetric tokamak plasma and the associated control systems [25,39,40]. The TSC simulates the time evolution of a free boundary plasma by solving the MHD equations on a rectangular computational grid. The MHD equations are coupled to the external circuits (representing poloidal field coils) through the boundary conditions. The code includes provisions for modeling the control system, external heating, and fusion heating.
The LSC module can also be called by the TRANSP [23,24] code. TRANSP represents the plasma with a axisymmetric, fixedboundary model and focuses on calculation of plasma transport to determine transport coefficients from data on power inputs and parameters reached.
This manual covers the basic material needed to use the LSC. If run in conjunction with TSC, the ``TSC Users Manual'' should be consulted [39,40]. If run in conjunction with TRANSP, online documentation will be helpful.
The few general mathematical subroutines (root finder, matrix invertor, random number generator) were taken from the the standard ``Numerical Recepies'' textbook [41]. X ray cross sections were entered from the literature [42,43]. Therefore no external libraries are needed for LSC.
A theoretical background of the governing equations and numerical methods is given in section 2. Information on obtaining, compiling and running the code is given in section 3.
This manual can be found at w3.pppl.gov through one or all of ~ignat, NTCC or topdac.
At such locations there will be links to the source code .
The manual is available in html written by Ian Hutchinson's program TtH and in html written by latex2html. For viewing TtH output with netscape under Unix you will need to put the following line in either .Xdefaults or .Xresources:
Netscape*documentFonts.charset*adobefontspecific: iso88591
Neither html version is meant for printing.
The manual is also available in PostScript as a .ps.gz file
and in the Portable Document Format pdf with activated internal links
as a .pdf file.
For LSC and TSC the poloidal flux is minimum in the center of the plasma and rises toward the plasma boundary. For LSC the B_{Z} at the outer midplane is positive for a positive current.
 (1) 
 (2) 
 (3) 
 (4) 
 (5) 
The antiHermitian part of the tensor is retained as a perturbation. For the case at hand, the principal such term enters as an imaginary correction to K_{zz}, and describes the interaction between the component of the wave electric field parallel to B and electrons whose speed along B matches that of the wave (Landau damping ). Thus, the plasma dielectric behavior is described by the following tensor elements, all other elements being zero:


 (13) 
In the coordinate system of equation (5), we decompose the dispersion relation (4) into its real and imaginary parts using equations (6)  (12),
 (14) 

If n º k c / w >> 1 a simplified dispersion relation can be found from the matrix equation leading up to equation (4) by asymptotically expanding

 (18) 
 (19) 
The extension of the local solutions to a spatially inhomogeneous plasma is accomplished by the eikonal method with the useful result that an initial wave field at r with an initial propagation vector k evolves according to Hamiltonian equations which preserve the local dispersion relation e_{r} = 0 along the ray trajectories

As in Hamiltonian mechanics, the spatial coordinates denoted by r are canonically conjugate to the wave number coordinates denoted by k. For example, in a Cartesian frame one would have (x,y,z) and (k_{x},k_{y},k_{z}) and fields varying as expi ò(k_{x} dx +k_{y} dy + k_{z} dz); and in a cylindrical frame one would have (R,Z,f) and (k_{R},k_{Z},l) where R, Z have dimensions of length, k_{R}, k_{Z} have dimensions of inverse length (wave number), and l is the dimensionless toroidal mode number. In the study of the axisymmetric tokamak, it is the cylindrical coordinate system that is the most natural since the toroidal wave number l is constant along the ray path ( d l / dt µ ¶e_{r} / ¶f = 0).
In formulating e_{r} it is necessary to use the canonical forms by writing k_{} = k ·B/B and k_{^}^{2} = k^{2}  k_{}^{2} where B and B are magnitude and vector quantities of the local magnetic field.
A spectral component of power W experiences a change in power DW over time interval Dt:

Note that equations (20), (21), and (22) also carry an index numbering the individual spectral component, which has been suppressed to improve readability. Given the velocity distribution and the profiles of macroscopic plasma parameters, the absorption of a lower hybrid spectrum can be computed. An actual incident wave spectrum is a continuous function of the parallel wave number. Computationally, this continuous spectrum is approximated by assigning the input power to a number of discrete rays as in equation (1), each ray having a definite initial k_{} and launched power. The number of rays and their distribution over k_{} can be chosen to adequately resolve the spectrum in cases of reasonably strong singlepass damping. On the other hand, if many passes are required to absorb the power, a very dense packing of rays into the spectrum might be required to obtain well behaved results.
The plasma is divided up into discrete shells between flux surfaces, typically 20  100. For each ray, equations (20) through (22) allow calculation of the needed quantities. The width of the shells and the detailed trajectory determine Dt in equation (22). As a practical matter, one cannot ensure that the shells are fine enough to keep DW / W < 1, so this quantity must be monitored and Dt adjusted downward the amount necessary to keep DW < W . The new Dt is recorded for use in equation (26) of section 2.3.1.
In applications of the model to many devices the following features are significant:
These effects are well established in the literature [14,18,19,20,21,32].
An important approximate relationship [18,19,21], exact under the electrostatic form of the cold plasma dispersion relation, between k_{0}, the the initial k_{} of a ray when launched, and k_{max}, the maximum possible value at a particular point in the plasma is the following.
 (23) 
Here R_{0} is the initial major radius of the LHCD ray of frequency w, and R, B_{p}, B_{f}, w_{pe} are the major radius, poloidal and toroidal magnetic field, and plasma frequency at some point of interest. The perpendicular dielectric constant is represented by S, a slowly varying quantity of order unity. It has been shown, without benefit of approximations [21], that this maximum tends to be reached in elongated plasmas after many reflections from the wall, whereas the maximum may not be reached in similar circular plasmas. It is worth noting that Kupfer, Moreau, and Litaudon [22] have used physics closely related to the basis of Eq. (23) to formulate a statistical theory of current drive with lower hybrid waves under conditions of multipass absorption.
 (24) 
 (25) 
 (26) 
 (27) 
 (28) 
We employ a onedimensional collision operator of the form,
 (29) 

In solving for f_{e} we set ¶/¶t = 0 because the time for equilibration between rf power and the electron distribution is short compared to the time for plasma to evolve. Then the solution for f_{e} is an integral in velocity space.
 (32) 
The rays do not move monotonically in ``radius'' (which in our calculation is measured by poloidal flux y ) and in fact typically make many inout changes of direction. To compute needed quantities at each flux surface, the motion of ray number i_{r} must be recorded by a zone number i_{z}, which is sequential with the intersection of that ray with any flux surface. For each zone index, the index i_{p} of the y surface intersected is kept, where i_{p} is 1 at the center of the plasma, increasing toward the edge. For each trajectory the i_{z} starts at 1 and increases to a maximum value; but i_{p} starts near its maximum, decreases initially, and exhibits individual behavior thereafter.
Changes in direction can be of several types. An inwardgoing ray can change to an outwardgoing ray, or viceversa, because:
speaking by analogy to a mechanical Hamiltonian system, the radial penetration limit is reached owing to the inpenetrability of a centrifugal barrier;
the root of the dispersion relation approaches zero (cutoff);
the ray is captured and specularly reflected near the plasma edge. This is accomplished by logic in the calculation which senses that a cutoff is approaching.
The quantities D_{ql;iv,ip} f_{e;iv,ip}, W_{iz,ir} are all made consistent by iteration, beginning with very small W_{iz = 1, ir} and ramping up toward the full power launched. We use an underrelaxation algorithm in the iteration of f_{e}, in that the past f_{e} is averaged with the newly computed f_{e} to get the new function propagated to the next cycle.
 (33) 
The key quantity is W_{s}(u), the energy (normalized to m_{e} v_{r}^{2}/2) imparted to the electric field E_{dc} by an electron as it slows down. The local electric field E_{dc} is either prescribed for a static simulation (usually set to zero) or supplied by TSC as part of an iteration to be described.
For very small electric field, W_{s}(u) can be represented in a power series. From the first two terms, one finds
 (34) 
The function W_{s}(u) behaves in such a way that for u approaching 1, corresponding to rf current drive cooperating with the ohmic current driven by the ambient electric field E_{dc} at a velocity near the runaway velocity, the calculated current grows extremely large. Because of the runaway condition, this is to be expected. Whenever u » 1, the rfdriven current is not properly calculable from our approach.
It is often necessary to iterate between TSC and LSC for the field and current at the exact time that LSC is called for power and current information, as follows. Suppose at time ^{(n1)} the local total current density is J^{(n1)}, the rfdriven current density is J_{rf}^{(n1)}, and the local electric field is E_{dc}^{(n1)}. We want to advance quantities to the new time, ^{(n)}, and be consistent with Ohm's Law: J = sE_{dc} + J_{rf}, where s is the parallel neoclassical conductivity. An initial guess at what the new electric field E_{dc}^{(n)} must be is formed from the new total current, as advanced by TSC, and the old rfdriven current and old electric field.
 (35) 
 (36) 
 (37) 
The new values E_{dc}^{(n)} and J_{rf}^{(n)} may not be consistent. This is checked by assigning them to the previous step labeled ^{(n1)} and repeating the process until the result is stable to a chosen accuracy. At that point, it is desired that the electric field is small enough that the runaway situation does not exist for electrons interacting with the waves. An inequality expressing this condition, as suggested by equation (20) of reference 12, is
 (38) 
Iteration is not necessary if when the ambient electric field is low and the density is high, such that w << v_{r} k_{}. In that case the evolution of the electric field and the current is handled without special attention as the transport calculation evolves.
Diffusing a velocity distribution is beyond the scope of the LSC, but we have introduced a heuristic approach similar to that used by V. Fuchs and others [17].
The idea is that the actual, diffused, RFdriven current J_{d} is modified from the computed RFdriven current J_{rf} by slowingdown, characterized by an inversetime n_{slow}, and a crossfield diffusion, characterized by D_{Jrf}.
Assuming for the moment slab geometry with x the spatial variable the diffusionlike equation we use is:
 (39) 
This can be transformed to a normalized flux coordinate [^x] º (y y_{min})/(y_{max}  y_{min}) with standard transformations for curvilinear coordinates. Under the reasonably good assumption that y ~ r^{2} the result is:
 (40) 
Equation 40 clearly has desirable properties: if D_{Jrf} is small then the J_{d} moves to J_{rf} on the n_{slow} time scale; and if D_{Jrf} is large then J_{d} and J_{rf} can be quite different, in particular J_{d} will be smooth compared to J_{rf} over a length of order Ö{ D_{Jrf} / n_{slow}}.
In LSC D_{Jrf} on input is one number, constant across the plasma, in units of meter^{2} per second. Internally, D_{Jrf} is an array, so generalization would be simple. The n_{slow} are found from equation 31 for a particular phase speed, corresponding to n_{} = 2, also constant across the plasma. The diffusion process is imagined to be complete at each call to LSC, so that the time derivative is set to zero. Under that condition, the solution for J_{d} comes from one inversion of a tridiagonal matrix.
The boundary condition on J_{d} at the outer boundary is a zero value. At the inner boundary, the center, the condition of zero flux (zero derivative with respect to radius) translates to the derivative with respect to [^x] being better behaved than 1/Ö{[^x]}, which cannot be expressed in a manner practical for the finite difference equations. Instead, we force dJ_{d}/d[^x] to be constant over the first two grid spacings. This is equivalent to specifying a zero value for a certain linear combination of J_{d} and dJ_{d}/d[^x] at the first interior grid point.
The rfdriven current found J_{rf} leads directly to the smoothed (diffused) J_{d} after the inversion of a tridiagonal matrix [41].
Note that the integral of J_{d} over the plasma cross section, or total diffused rfdriven current, is not constrained to have any particular relationship with the same integral over J_{rf}, the total undiffused current. In fact, it seems reasonable that diffusion can in some circumstances increase current by moving fast electrons to regions of lower collision rate. However, it appears that the diffused current is less than the undiffused current for many actual situations, typically having deposition at midradius and mildly peaked density profile.
The heuristic smoothing method of V. Fuchs [17] is similar to the method described here, but apparently treats D_{Jrf} and n_{slow} as constants, sets the change in total current from an estimate incorporating the formula for n_{slow}, and finds the value needed for D_{Jrf} as an eigenvalue.
Our model is not a simulation of energetic electron diffusion in that it does not transport fast particles or change the rf wave damping because of their motion.
The undiffused current J_{rf} was computed from the ray tracing, quasilinear development of an electron distribution function in the parallel velocity giving the local power deposition per unit volume P_{0}, and finally the KarneyFisch [12] electric field adjustment to current density.
As explained in the previous section, the diffused current density J_{d} can be quite different from J_{rf}, and therefore quite different in radial distribution from the P_{0} found from ray tracing. This situation is contrasts with the intuition that power density is proportional to current density times background number density. Therefore, we added an option to spread the deposited rf power according to
 (41) 
The PBXM experiment to which we have applied our calculation the most has a 2dimensional xray camera with an axis roughly tangent to a midplane circle just inside the major radius of the the plasma. The images found by this camera reveal clues to the behavior of the fast electrons in the plasma. At the same time, such images can be computed from the electron velocity distribution we establish in our model.
Taking account of obstructions, the actual location of pinhole, focal plane, etc., we evaluate the 2d image with the following formula.

In the above,
I_{xp}
is the xray intensity on pixel p of area dA_{p}
at the camera focal plane,
dl is an element of path length
in the plasma on a line passing through
the pinhole and the pixel,
dA is an area element for emission into the pixel p,
dW_{p} is the solid angle subtended by the pinhole
as seen from the emitting
volume element,
n_{i(e)} is the ion (electron) density at the emitting location,
E¢ denotes photon energy emitted by an electron of energy E,
¶^{2} s_{ei}/¶E¢¶W is the cross section for
electronion bremsstrahlung differential in energy and solid angle,
[^(n)] is a unit vector connecting the emitting volume with
the camera element, and
[^(b)] is a unit vector along the magnetic field, and
T(E¢) is the transmission factor through vessel windows
and other absorbers.
The differential cross section used in our calculation is that given
by formula 2BN of Koch and Motz
[42].
The significance of ``2'' is the differential in energy and angle;
``B'' is Born approximation;
``N'' is no screening of the ions.
Note that the distribution function f_{e} is that computed from
the onedimensional FokkerPlanck model, and is itself not influenced by the
electric field E_{dc}.
The absorption from vessel and windows is computed from fits
to experimental absorption data
[43].
There are graphics calls embedded in the code for PPPL's ``Scientific Graphics,'' or SGlib , or just SG, and for ``gnuplot.''
Information and source for SGlib is at:
w3.pppl.gov/rib/repositories/NTCC/catalog/Asset/sglib.html.
Information and source for gnuplot is at:
www.cs.dartmouth.edu/gnuplot_info.html.
The combination LSC/gnuplot should work on any Unix, VMS, MSDOS, or Windows system with a Fortran 77 compiler. Gnuplot for he Macintosh has been reported, but not found on the Internet in September 2000.
The combination LSC/sglib should work on any Unix or VMS system with a Fortran 90 compiler and a C compiler.
The source for LSC is constructed with calls for both sglib and gnuplot imbedded, but there are exceptions, as follows. Contour plots and the 3Dlike isometric sketch of ray paths in the toroid require sglib to be loaded; and, none of the graphs used by the X ray camera have been converted to gnuplot. Converting contour plots to gnuplot does not appear to be possible in a satisfactory way with a modest amount of work, and will probably never be done. Fortunately, these plots never seemed to be particularly useful. Converting the 2D plots of the X ray camera to gnuplot is probably easy enough, and can be done if there is a demand for these graphs under gnuplot.
By default all calls to sglib are compiled into the executable code as nooperations with the cpp command #ifndef USE_SGLIB. In a similar way the X ray camera is not compied via the cpp command #ifndef USE_X_RAY_CAMERA.
As mentioned in the introduction, this manual and links to the source code can be found at w3.pppl.gov, through one of ~ignat, NTCC or topdac . The source code should be accompanied by at least one ``Makefile'' which works in the development environment, Linux Fujitsu Fortran, for the standalone code. Linux GNU Fortran 77 (g77) was also used in the development, however, the g77 executable is somewhat slower. LSC is constructed as a subroutine to be called as a module of a transport code such as TSC or TRANSP, but standalone use with a simple driver program is easily done. The distribution provides a driver, lscdrive.F , and a simplifed mockequilibrium writer, WrCircEq.F . Standalone runs can be extremely useful in studying the behavior and assuring reasonable convergence with the parameters chosen.
Naturally, standalone computations cannot show the interaction of the current driven with the electric field generated, reacting to the current, in a complete calculation.
Details follow in the following subsections.
A header line (informational comments) equhd is read first followed by integers npsitm, nspc, kcycle, a plasmatime in the calling program, times (which is ignored by LSC) plasma parameters anecc, tekev, anicc, tikev, amass, achrg, the electric field (loop voltage) voltlp, and specifications of the equilibrium, rho, vptemp, xsv, pary, ppary, gary, gpary, nx, nz, isym, iplim, psimin, psilim, xmag, zmag, rgzero, bgzero, apl.
The meaning of the variables, and units used, are documented in the source.
The specification is done through the NAMELIST (listdirected) file input.lhh under the namelist name inpvalue
The namelist variables are documented in source file Inpval.inc .
The number of rays to be launched is nrays, which decomposed into the product of ntors and npols, so each toroidal spectral component can have a number of poloidal components, each of which behave somewhat differntly. Defaults are in params.inc .
There are two methods of representing the spectrum, a simple model of Gaussian shapes, and a calculation stemming from Brambilla paper [1].
Three tables list the LH parameters according to the following plan.
Table 1 lists the most important input paramters, and also the ones needed to specify a Gaussian spectrum.
Table 2 lists the input parameters to choose a Brambilla calculation of the waveguide spectrum. Note that some of the parameters of Table 1 are ``recycled'' for the Brambilla case, and therefore have a somewhat different meaning: nGrps, powers.
Table 3 lists the input parameters for the heuristic model of current and power diffusion.
fghz  The launched LH wave frequency (in GHz), typically 2.45 or 4.6 (GHz). 
nGrps  The number of groups (or couplers, which see) used in calculating the launch spectrum; the default is NGRPDIM (3). 
centers()  The location of the peaks in the Gaussianmodel spectrum, in n_{}; the three defaults are all 4.0 . 
widths()  The widths of the model spectrum peaks, as Dn_{}; the three defaults are all 1.0 . 
powers()  The relative powers in either the Gaussianmodel spectrum components (if DoBram is zero) or in the couplers (if DoBram is equal to unity); as an example, if powers(1) and powers(2) are both 1.0, LSC will assign equal powers to the first two components/couplers, whereas if powers(1) = 1.0 and powers(2) = 0.5, LSC would assign twice as much total power to the first component/coupler as the second. The default values are 1.0, 0.1 and 0.1 for elements 1, 2 and 3, respectively. 
nrays  The number of rays LSC will trace; its most natural value (also the default) is the product of ntors with npols, but if nrays is given as something different, then a sensible action is taken. 
ntors  The number of distinct toroidal components in the number of rays (nrays); the default setting for ntors is NTORDIM (30). 
npols  The number of poloidal variations to each toroidal variation; in part, this is a way to increase the number of rays launched from the same Bramilla spectrum, but also effective for the model spectrum. The default 1. 
nparmax, nparmin  The floating point numbers giving the maximum and minimum values of n_{} over which to distribute the ntors number of toroidal values when using the Maxwellianmodel spectrum (hence, these inputs are ignored when DoBram is set to 1). The defaults are 5.5 for nparmax and 2.5 for nparmin. 
npolmax, npolmin  The floating point numbers giving the maximum and minimum values of n_{q} over which to distribute the npols number of poloidal values when using either the Brambilla or Gaussian spectrum; their default values are 1.0 for npolmax and 1.0 for npolmin. 
DoBram  INTEGER: If DoBram is equal to 1 (the default), then LSC will perform the Brambilla calculation, using nGrps as the number of couplers, and the array powers() as the relative distribution of power between the couplers. If DoBram is 0, LSC uses Gaussian models for the launched spectrum. 
couplers()  The dimensions of certain couplers are built into LSC in DATA statements; hence, all the user must specify is the eightcharacter name(s). Presently, LSC will recognize PBXMFAST, PBXMSLOW, TOKDEVAR, SLOWSLOW, TORSUPRA, TFTRLHCD, JET_LHCD, USRSPEC1, USRSPEC2, USRSPEC3 . 
phaseDeg()  The relative phases (in degrees) applied to the couplers; used in computing the Brambilla spectrum; the default values are 90.0, 180.0 and 135.0, respectively. 
nslices  The number of slices on the Brambrilla spectrum. The default is 301. 
Regarding Table 2, the `USRSPECn' waveguides are set equal to TFTRLHCD in the code as distributed. The intent is for the data statements to be replaced as needed. The format of the specification of the waveguide coupler is the following: az(i) is the toroidal distance to the near edge of i^{th} guide in cm; bz(i) is the toroidal distance to the far edge of i^{th} guide in cm; so that az(2)bz(1) is the the width of the septum between the first two guides.
DiffuJrf  The heuristic diffusivity for rfdriven current in meter squared per second. It balances a collision rate, as shown in equation 40. The default is 0. 
PrfSpred  This governs the heuristic spreading of rf power as given by a in equation 41. Dimensionless, between 0.0 and 1.0. The default is 0. 
Table 4 gives the computational parameters, and lists the defaults for each one.
HstpLH  The length of the ray integration step in meters. The default is 0.005 (meters). 
nstep  The maximum number of integration steps (each having length HstpLH) through which LSC will track a ray; the default is 20,000 steps. 
npsi  The number of y bins used in LSC (during the D_{ql} calculation); npsi can be either more or less than the npsii from the TSC grid. The default is NPSIDIM (100). 
nzones  The number of zone (psi surface) crossings through which LSC will trace a ray. The default is NZONDIM (2000). 
nv  The number of parallel velocity bins; the default is NVELDIM (401). 
nsmoo  The number of velocity bins over which smoothing is to be done. The default is NSMODEF (9). 
nsmw  The characteristic width (in number of velocity bins) of the smoothing. The default is NSMWDEF (3). 
nRampUp  The number of steps taken in ramping up the power for the iteration of f_{e}, D_{ql} and P_{rf}; some of these steps may be with a constant power (controlled by nFlat, above). The default is NRAMPDIM (200). 
nFlat  The number of rampup iterations with a flat (constant) power. The default is NFLATDEF (10). 
WeghtItr  The proportionality of the weighting of the old and new electron distribution functions during each rampup iteration. The default is 0.20  20 percent new f_{e} and 80 percent old f_{e}. 
Do1Rpr  This INTEGER flag modifies the way LSC recalculates rays when it is being called from TSC or TRANSP. If Do1Rpr is equal to 1, then LSC recomputes only 1 ray per call from TSC/TRANSP after the initial computation of all rays; this smooths the dynamic behavior. If equal to 0, then LSC recomputes all rays when called by the driving program. 
ScatKdeg  The angle in degrees by which the perpendicular component of the wave vector, k_{^}, may be rotated on rms average when the ray ``bounces'' (i.e. as in off of the wall). If ScatKdeg is 0.00 (the default), then there is no scattering  the ray is simply reflected specularly. 
TurnNegs  If this INTEGER flag is set to one, then any ``negativelydirected'' Brambrilla spectrum elements will be turned around; literally, the corresponding n_{} components will be multiplied by a negative one, which will reverse the direction of the spectrum element. The default is zero. 
nfreq  After each ray has been traced for an nfreq number of steps, LSC takes stock of what has occurred. This INTEGER has a default value of 100 (steps). 
PlFlg()  This INTEGER array is made up of the flags that inform LSC as to which (if any) of the possible plots should be produced as output; if a particular plot is desired, then the corresponding PlFlg element should be set to one (1). The default settings are all zeroes (no plots). For a full description of the possible output plots, see the section below on LSC Outputs. 
PrFlg()  This INTEGER array is similar to PlFlg, above, but holds the flags governing LSC ``screen'' output. LSC has three output channels which have their destinations by the driving program, and conceivably any or all of these channels could be linked with the screen. Alternatively, these channels could be linked to additional output files. For a more complete description of the channels and the possible outputs, see the section below on LSC Outputs. The default setting of this array is all zeroes. 
idiag()  This INTEGER array specifies at which rampup iterations the results of f_{e} and D_{ql} calculations will be ouputted; that is, the first element of idiag tells LSC at which iteration to produce the first set of plots, the second element tells at which iteration to produce the second set, etc. idiag starts out with all IDIAGDIM elements set to zero; the user should (at the very least) set the first element of idiag to nRampUp  otherwise, LSC does not know ``when'' to produce certain outputs, and will not (regardless of the values set in PlFlg). 
Typical input namelists, to be found in input.lhh, are found in Table 7.
&inpvalue fghz = 4.6, nGrps = 1, powers(1) = 1.00, powers(2) = 1.00, phaseDeg(1) = 130., phaseDeg(2) = 130., nslices=301, ntors=15, npols=1, DoBram = 1, couplers(1)='TFTRLHCD', couplers(2)='TFTRLHCD', DiffuJrf = 0.00, / &inpexprt HstpLH = .01, nstep = 20000, nfreq = 50, npsi = 100, nzones = 2000, nv = 401, nsmoo = 9, nsmw = 3, weghtitr = .2, plflg( 1) = 0, 1, 0, plflg( 4) = 1, 0, 0, 0, prflg( 1) = 0, 0, 0, incLabls = 1, /
Table 8 lists the arguments that must be supplied by the calling program, and is largely selfexplanatory. The power is the most significant parameter. For use of LSC with a transport code and determination of both current and electric field in a selfconsistent way iRayTrsi is crucial, and almost always equal to 0 because of the number of iterations needed.
Wr2TSC.out  Electron power, electron current, E/J dJ/dE, ion power 
LSCcomm.out  
LSCcom2.out  
ray.dat  
Name of Module  Brief Description of Role in the Computation 
AcDc.F  Main; changelog; block data; some i/o 
Rayini.F  Ray initialization 
Rayio.F  Ray i/o 
Raystore.F  Ray information stored 
Raytrace.F  Ray tracing 
Wr2TSC.F  Write power and current for calling program; other output 
WrCircEq.F  Write out a `practice' circular equilibrium 
XrayCam2.F  Predict results for an X ray camera 
brambJES.F  Compute a launced spectrum from the Brambilla approach 
cycle.F  Ramp up power to get electron distribution 
ezsg.F  Scientific Graphics (SG) shell subroutines 
fe.F  Initialize f_{e}; do smoothing and derivative 
fits.F  Fits to current response from the KarneyFisch approach 
gnup.F  Gnuplot shell subroutines 
grap.F  Grid interpolations 
grapgrd.F  Find quantities by interpolation 
gridgen.F  Generate grids 
grids.F  Generate grids, also 
jrf.F  Form rf current 
lscdrive.F  Call LSC in standalone mode 
matr.F  Manipulate matrices 
pbxio.F  Output, mostly from PBX days 
plasejv.F  Find plasama parameters (EJ Valeo) 
power.F  Develop power deposited by rays 
ql.F  Develop D_{ql} 
Name of Module  Brief Description of Role in the Computation 
CGSetc.inc  Constants in the cgs system 
Doflags.inc  Flags for dooptions, such as DoBram 
DqlBins.inc  Quantities for D_{ql} 
Escan.inc  Userinput E_{dc} for standalone testing 
FeBins.inc  Quantities for f_{e} 
Implic.inc  `IMPLICIT NONE' and #define s for CPP 
Inpval.inc  Namelist inputs 
Jrf.inc  Quantities for J_{rf} 
MKSetc.inc  Constants in the mks system 
PIetc.inc  Constants related to p and roots 3/2, 1/2 
PlPr.inc  Flags for PLotting and PRinting 
ProfBody.inc  Quantities specifying the plasma 
Ramppwr.inc  Parameters and quantities governing rampup of power 
RayBins.inc  Quantities for storing ray information 
RayWrk.inc  Quantities govening ray calculation 
TSCgrap.inc  Quantities specifying the plasma 
WkAry.inc  Scratch and miscellaneous arrays 
dielec.inc  Dielectric tensor elements; quantities for the present ray 
gnuI.inc  Integers defining gnuplot graphs 
gnuR.inc  Reals defining gnuplot graphs 
params.inc  Parameters specifying array sizes 
plcmx.inc  Steps in R and Z for the equilibrium grid 
power.inc  Quatities for power deposited 
sgIbk.inc  Integers defining sglib graphs 
sgRbk.inc  Reals defining sglib graphs 
tscunits.inc  Unit numbers and flags from calling program 
Name of Module  Brief Description of Role in the Computation 
Xray.inc  Quantities describing xray properties: foils, energy bins, etc 
camera.inc  Quantities describing the camera 
emitter.inc  Quantities describing the plasma emitter 
emparams.inc  Quantities describing the emission calculation 
numerics.inc  Quantities in common for the xray camera 
xgraphs.inc  Working arrays 
xparams.inc  Parameters for computational bins 
&inpxry E_max = 0.20, E_min = 0.01, E_ph_min = 1.0e6, FoilCode = 'AG', PusherMajor = 1.152, PusherMinor = 0.190 R_bound_max = 2.0, R_bound_min = 1.1, R_tangent = 1.525, z_tangent = 0.0, Rpinhole = 2.661, Zpinhole = 0.0, pinhole_size = 0.00625, pinhole_size = 0.00625, dE_ph = 0.01, dFoilTCM = 0.0, dlxray = 0.005, focal_length = 0.495, iAbsXray = 1, screen_d = 0.215, /
E_max/min  To render the calculations tractible, only the contributions from electrons in the plasma with kinetic energies within a certain range will be considered; E_min and E_max (REAL numbers) specify (in keV) the minimum and maximum electron energies of this range. If the user does not specify these energy bounds, LSC will set E_min to 0.01 keV and E_max to 2.0 keV. 
E_ph_min  E_ph_min specifies the minimum energy that a photon must have in order for LSC to add its contribution into the line integral. The default value of E_ph_min is 1×10^{6} keV. 
FoilCode  FoilCode is a twoCHARACTER parameter specifying the type of metal foil that the XRays must pass through in going from the pinhole to the film; FoilCode can be one of the following: ``CU'', ``AG'', ``TA'', ``MO'', or ``00'' for no foil at all. By default, FoilCode is set to ``AG''  a silver foil. 
PusherMajor  PusherMajor, PusherMinor, are the radial positions of the major and minor ``pusher'' field coils, as measured from the axis of toroidal symmetry. Their default values are 1.152 meters and 0.190 meters, respectively. 
R_bound_max  R_bound_max, R_bound_min, specify the radial extent of the plasma; that is, the plasma will extend no further in than R_bound_min, and no further out than R_bound_max. The default values for R_bound_min and R_bound_min are 1.10 meters and 2.00 meters, respectively. 
Z_bound_max  Z_bound_max, Z_bound_min, are like R_bound inputs above, Z_bound_min and Z_bound_max specify the furthest extents of the plasma along z. If the user does not specify them, Z_bound_max and Z_bound_min are set to 1.00 meters and 1.00 meters. Basically, these six inputs specify the boundaries at which LSC would stop calculating the line integrals, since there is no plasma beyond the R or Z _bounds or next to the ``pusher'' coils. 
R/z_tangent  The pinhole aperture of the XRay camera and the point at the center of the image screen define a line in space; clearly, this line will be tangent to some circle about the axis of toroidal symmetry; R_tangent is the radius of this circle (with a default value of 1.524 meters), while z_tangent is the zposition (height) of this circle (with a default value of 0.0 meters). 
R/zpinhole  These three inputs completely specify the coordinate location of the pinhole within the torus. The default values for these three are 2.661 meters, 0.0 meters, and 0.0 degrees, respectively. 
dE_ph  The interval between photon energy bins, in keV. The default is 0.01 keV. 
dFoilThickCM  The thickness of the metal foil, in units of centimeters. The default value is 0.00 centimeters. 
dlxray  The default value of this REAL parameter is 0.005 meters. 
focal_length  This is the focal length in meters of the XRay camera  literally, this is the distance from the pinhole to the screen upon which the image is formed. The default value is 0.4905 meters. 
iAbsXrays  This INTEGER flag tells LSC whether or not to take account of the absorption of Xrays by the foil when calculating the Xray camera image; the default value of iAbsXrays is 1 (absorption). 
nEbins  This INTEGER specifies the number of XRay energy bins. If left unspecified, the default value is NENDIM. 
nMUbins  This INTEGER specifies the number of mu bins used in the calculation of the line integral through the plasma to the pinhole. By default, nMUbins is set to NMUDIM. 
n_pixel_x/y  These REAL input parameters are (respectively) the extent of the XRay camera image, in the x and ydirections, in number of pixels; that is, the camera image will be n_pixel_x pixels "long" along x, and n_pixels_y "long" along y. Their default values are both NPIXDIM. 
npoints_max  By default, npoints_max is set to MAXPOINTS. 
nr/z_source  Are the number of sections the plasma crosssection is divided into, radially and heightwise. The default settings for nr_source and nz_source are NRDIM and NZDIM, respectively. 
pinhole_size  The physical size (diameter in meters) of the pinhole; the default is 0.00625 meters. 
screen_d  The "diameter" in meters of the screen upon which the XRay camera image is formed; screen_d is the diameter of the largest circle that can be inscribed in the screen. By default, screen_d is set to 0.215 meters. 
F. W. Perkins devised the approach taken in the consideration of scattering. M. Ono , F. Paoletti , and H. Takahashi contributed discussions regarding the limited upshift of parallel index of refraction. S. von Goeler , S. E. Jones , and D. P. Enright provided valuable help regarding the xray camera. V. Fuchs encouraged the approximate treatment of rf current diffusion.
Operation of the code by M. Shoucri (of the Tokamak de Varennes) in conjunction with TSC and by S. M. Kaye , R. V. Budny (of PPPL) and P. Stubberfield (of JET) in conjunction with TRANSP led to early improvements.
Operation by M. Ju (of KSTAR) led to a recent and significant improvement in the stability of the contruction of the electron distribution function from ray data. S. Bernabei and M. Okabayashi supported, encouraged, and advised this project over the past several years.
An earlier version of this manual depended on the work of
A. J. Redd
. The present manual reuses
many of his contributions.