Lower hybrid Simulation Code (LSC) Manual by DW Ignat, SC Jardin et al

Lower Hybrid Simulation Code Manual

# Lower Hybrid Simulation Code Manual

## Abstract

The Lower hybrid Simulation Code (LSC) is computational model of lower hybrid current drive in the presence of an electric field. Details of geometry and the plasma profiles are treated. Two-dimensional velocity space effects are approximated in a one-dimensional Fokker-Planck treatment. If LSC is coupled with a model that connects local toroidal electric field and local current through Faraday's law, then details of current development can be examined. This Fortran 77 code does not call mathematical subroutines from proprietary libraries, and has been run under several compilers on several platforms. Graphical output is available via the public domain Scientific Graphics library developed at PPPL (sglib) and the gnuplot system, which is also freely available. LSC is therefore portable. The intent is to operate in double precision so results are not sensitive to changes between machines with a natural 64 bit word and 32 bit work stations, but operation in single precision is possible. Operation in double precision does not depend on compiler switches (for example, -r8) but is set up for the use of #define REAL real*8 and the standard C-pre-processor.

To read the mathematics in this version, which was produced by TtH, under Unix and Netscape, put the following line in .Xdefaults or .Xresources:

Nov 19, 2000

# Contents

1  Introduction
2  Theoretical Background
2.1  Geometrical Configuration
2.2  Wave Propagation
2.2.1  Dispersion relation and energy absorption
2.2.2  Changes to k||
2.3  Fokker-Planck Equation
2.3.1  Quasilinear Diffusion Coefficient
2.3.2  Kinetic Equation
2.3.3  Some Numerical Details
2.4  RF-Driven Current
2.6  Heuristic Power Diffusion Estimate
2.7  X-ray Camera Image
3  LSC Usage
3.1  Input Files
3.1.1  Plasma
3.1.2  LH parameters
3.1.3  Computational parameters
3.1.4  Flags for graphs and print files
3.1.5  Calling arguments
3.2  Output files
3.3  Output graphs
4  Appendix: Brief Description of Source Files

## 1  Introduction

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]. Two-dimensional velocity space effects are approximated in a one-dimensional Fokker-Planck 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, fixed-boundary 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, on-line 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:

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.

## 2  Theoretical Background

The following description of the governing equations and numerical methods closely follows the original paper describing the LSC model [32].

### 2.1  Geometrical Configuration

A cylindrical coordinate system with R as the radial coordinate, Z as the axial coordinate and f as the toroidal symmetry angle is used. In LSC this triple is thought of in as (R,Z,f), whereas in TSC this is considered (R,f,Z). In other words, there can be confusion over the f direction.

For LSC and TSC the poloidal flux is minimum in the center of the plasma and rises toward the plasma boundary. For LSC the BZ at the outer midplane is positive for a positive current.

### 2.2  Wave Propagation

#### 2.2.1  Dispersion relation and energy absorption

We assume that the plasma varies sufficiently slowly on a wave period and a wave length so that the WKB approximation is valid. We therefore assume that the wave electric field can be decomposed into a set of components indexed by j,

 E = Sj Ej(r ) expi Fj  ,
(1)
where the rapid space and time variation occurs through the exponential

 Fj º óõ kj(r ) ·dr  -  wt  ,
(2)
which depends on the local wave vector kj and the frequency w. Additional time behavior is ignored because the transit time of a wave is short compared to the time scale of plasma evolution. The mode amplitudes Ej vary much less rapidly. Suppressing the index j, each mode locally satisfies the matrix equation, [44]

 [ kk - I k2 + k02 K(r,k, w) ]·E = 0
(3)
where k0 = w/c represents the free-space wavenumber, and K is the dielectric tensor, which has a nontrivial solution for a vanishing determinant

 e(w, k, r) º |[ kk - I k2 + k02 K(r,k, w) ]| = 0     .
(4)
The dispersion relation  (4) connects the vector k and frequency w at the spatial location r (explicit time variation of the dielectric properties are assumed negligible from here on). To proceed it is convenient to choose a local Cartesian coordinate system [44], such that [^(z)] ×B = 0, and k is contained in the x-z plane:

 k = k|| ^z + k^ ^x .
(5)
For lower hybrid waves, the ion cyclotron frequency is much less than the wave frequency which is much less than the electron cyclotron frequency: wci << w << wce, an inequality assumed throughout. In the Hermitian part of the dielectric tensor we keep only cold plasma terms, except that the dominant warm-plasma term is carried to guard against singularity near the lower hybrid resonance.

The anti-Hermitian part of the tensor is retained as a perturbation. For the case at hand, the principal such term enters as an imaginary correction to Kzz, 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:

 Kxx = Kyy
 =
 S - ak^2     ,
(6)
 -Kyx = Kxy = i D
 =
 i wpe2wwce ,
(7)
 Kzz
 =
 P     + i Kzz,i    ,
(8)
where

 S
 =
 1 + wpe2wce2 - å j wpi,j2w2 ,
(9)
 a
 =
 34 wpe2wce4 vTe2 + 3 å j wpi,j2 vTi,j2w4 ,
(10)
 P
 =
 1 - wpe2w2 ,
(11)
 Kzz,i
 =
 - p wpe2w óõ dv||v|| ¶fe¶v|| d(w- k||v||)     ,
(12)
and the new symbols are as follows: wpe is the electron plasma frequency ne e2/(e0 me), wpi,j is the ion plasma frequency of the jth species, vTe and vTi,j are electron and ion thermal velocities ( º kT/m), and fe(v||) is a one-dimensional electron velocity distribution function normalized such that

 óõ dv||fe(v||) = 1  .
(13)
The permittivity of free space is e0. The thermal term designated by a would be important near lower hybrid resonance , S = 0. However, in cases treated in this paper S does not approach zero, and a is not important. This would be true whenever w2 >> wcewci. The wave-particle interaction responsible for electron heating and current drive is in Kzz,i. In the event the lower hybrid resonance should become important, a thermal term representing the ion wave-particle interaction would have to be added to Kxx and Kyy, but here we assume such terms vanish.

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

 e = er + i ei = 0  ,
(14)
where

 er
 =
 -ak^6 + k^4 S + k^2 éë (P + S)(k||2-k02S) +k02D2 ùû
 +P éë (k||2-k02S)2-k04D2 ùû ,
(15)
 ei
 =
 ¶er¶P Kzz,i     .
(16)

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

 ^n
 ~
 ^n 0 + 1n2 ^n 1 + O( 1n4 )  ,
 E
 ~
 E0 + 1n2 E1 + O( 1n4 )  .
(17)
To lowest order, I^·E0 = 0, where I^ = I - [^(n)] [^(n)], so that E0 = [^(n)]0 E0. In first order, there arises the solubility condition

 ^n 0 ·K · ^n 0 = 0  ,
(18)
which gives the electrostatic limit of the dispersion relation :

 er = - ak^4 + k^2 S + k||2 P     .
(19)
The consistency condition |E1| = |K ·E0| < n-2 |E0| is satisfied for n||2 > 1 + (wpe2 / wce2).

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 er = 0 along the ray trajectories

 drdt
 =
 - ¶er¶k / ¶er¶w
(20)
 dkdt
 =
 + ¶er¶r / ¶er¶w
(21)

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 (kx,ky,kz) and fields varying as expi ò(kx dx +ky dy + kz dz); and in a cylindrical frame one would have (R,Z,f) and (kR,kZ,l) where R, Z have dimensions of length, kR, kZ 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 µ er / f = 0).

In formulating er it is necessary to use the canonical forms by writing k|| = k ·B/B and k^2 = k2 - 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:

 DW
 =
 -2  ei / ( ¶er¶w ) W Dt
 =
 -2 ¶er¶P Kzz,i /( ¶er¶w )  W Dt
 =
 2p wpe2w óõ dv||v|| ¶fe¶v|| d(w- k||v||) ¶er¶P /( ¶er¶w )   W Dt     .
(22)

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 single-pass 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.

#### 2.2.2  Changes to k||

In applications of the model to many devices the following features are significant:

1. toroidal effects move the value of k|| up and down during propagation - a rise helping absorption, and a fall retarding absorption
2. the rise in k|| is a stronger effect as density, poloidal field inside the plasma increase, and as aspect ratio (R/a) decreases

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.

 k||max £ k||0 × R0/R1-(Bp/Bf)( wpe/w) S-1/2 .
(23)

Here R0 is the initial major radius of the LHCD ray of frequency w, and R, Bp, Bf, wpe 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 multi-pass absorption.

### 2.3  Fokker-Planck Equation

#### 2.3.1  Quasilinear Diffusion Coefficient

An incremental contribution to the quasilinear diffusion coefficient Dql at velocity v|| from a wave field of wave number k|| is given by

 Dql(v||) = p2 æç è eme ö÷ ø 2 E||2 d(w- k|| v||),
(24)
where E|| represents the amplitude of the wave field parallel to the local static magnetic field, and the electron charge and mass are e and me [44]. One instructive way to find the relationship between field E|| and wave power W is to equate Pql, the energy per unit time per unit volume going into electrons and out of the wave from the quasilinear point of view,

 Pql = - ne me óõ dv||v||Dql (v||) ¶fe¶v||
(25)
with the similar quantity from the ray point of view as assembled from equations (12), (16), (22) to obtain

 Dql(v||) = 2 æç è pe0 ö÷ ø æç è eme ö÷ ø 2 W ¶er /¶Pw¶er/¶w æç è DtDV ö÷ ø d(w- k|| v||)
(26)
for the incremental Dql from a wave of power W traversing a flux shell of volume DV in time Dt.

#### 2.3.2  Kinetic Equation

An electron kinetic equation can be written

 ¶fe¶t = æç è ¶fe¶t ö÷ ø c + æç è ¶fe¶t ö÷ ø w
(27)
Here the ()c term is the Coulomb collision operator, and the ()w term is the wave diffusion (quasilinear) operator. There is no term involving a steady electric field. The wave diffusion operator is the one-dimensional divergence of the rf-induced flux:

 æç è ¶fe¶t ö÷ ø w = ¶¶v|| Dql(v||) ¶fe¶v|| ,
(28)
where now the term Dql signifies a sum over all waves in existence on a flux surface, with the appropriate powers and velocities. Note that the use of a simple sum means that we assume there are no interference effects.

We employ a one-dimensional collision operator of the form,

 æç è ¶fe¶t ö÷ ø c = ¶¶v|| éê ë æç è Dc(v||) ¶¶v|| + nc(v||)v|| ö÷ ø fe(v||) ùú û .
(29)
The collisional diffusion and drag coefficients are given by

 Dc(v||)
 =
 bZ  G/vTe  [1 + (v||/ vTe)2]-3/2     ,
(30)
 nc(v||)
 =
 bZ  G/vTe3  [1 + (v||/ vTe)2]-3/2     ,
(31)
where G = lnL  ne e4 / (4pe02 me2 ) and the normalization coefficient bZ is chosen to yield the correct value for the electrical conductivity in the absence of rf. Approximately, bZ = (1 + Z) / 5, with Z the effective ion charge. The collision operator has several correct properties, such as conservation of particles, producing a Maxwellian (with thermal velocity vTe) in the absence of wave power, and correct asymptotic velocity dependence of the coefficients for high speeds v >> vTe.

In solving for fe 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 fe is an integral in velocity space.

fe(v||) = 1
 Ö 2pvTe2
exp æ
ç
è
- ó
õ
v||

0
nc(v¢)  v¢ dv¢
Dc(v¢) + Dql(v¢)
ö
÷
ø
(32)

#### 2.3.3  Some Numerical Details

Our calculation requires a discrete grid in v|| indexed by iv and many discrete rays indexed by ir which at each flux shell crossing indexed by iz have a definite w/ k|| and a power W at a definite location in the plasma indexed by ip. The delta function d(w- k||v||) connects the discrete velocity quantity with the discrete wave number quantity and therefore must be interpreted to have a width that causes smooth variations in velocity while the spectrum of waves, represented by a limited number of rays, is not smooth. The width of the delta function in velocity should be such that for two nearby rays separated by Dk|| the range in velocity Dv|| is approximately the relative separation of two nearby rays Dk||/k|| times the phase velocity w/ k||. However, some judgment should enter the actual selection of the width of the smoothing function because k|| and therefore Dk|| can evolve considerably. If the smoothing range is too narrow in velocity, then non-physical holes'' can appear in the Dql with a resulting ragged behavior of fe. On the other hand, if the smoothing range is too great, then a wave of phase speed w/k|| interacts with electrons of much different parallel velocity v||, especially much lower parallel speed, since there are many more low speed electrons than high speed electrons. This can lead to spurious power absorption.

The rays do not move monotonically in radius'' (which in our calculation is measured by poloidal flux y ) and in fact typically make many in-out changes of direction. To compute needed quantities at each flux surface, the motion of ray number ir must be recorded by a zone number iz, which is sequential with the intersection of that ray with any flux surface. For each zone index, the index ip of the y surface intersected is kept, where ip is 1 at the center of the plasma, increasing toward the edge. For each trajectory the iz starts at 1 and increases to a maximum value; but ip starts near its maximum, decreases initially, and exhibits individual behavior thereafter.

Changes in direction can be of several types. An inward-going ray can change to an outward-going ray, or vice-versa, because:

1. the roots of the dispersion relation converge (mode conversion);
2. speaking by analogy to a mechanical Hamiltonian system, the radial penetration limit is reached owing to the inpenetrability of a centrifugal barrier;

3. the root of the dispersion relation approaches zero (cutoff);

4. 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 Dql;iv,ip fe;iv,ip, Wiz,ir are all made consistent by iteration, beginning with very small Wiz = 1, ir and ramping up toward the full power launched. We use an under-relaxation algorithm in the iteration of fe, in that the past fe is averaged with the newly computed fe to get the new function propagated to the next cycle.

### 2.4  RF-Driven Current

We calculate the current driven on each flux surface according to our equation (33) which follows the prescription given in equation (21b) of reference (12), except that we drop the term arising from a non-zero runaway probability , obtaining

 Jrf = -e nenr óõ dv||Dql(v||) ¶fe(v||)¶v|| · ¶Ws(u)¶u .
(33)
In the above, nr = G/|vr|3, u = v||/vr, and vr = - sign(eEdc) Ö{meG/|eEdc|}, as explained in the reference [12]. These definitions give as Jrf the current density of stopped electrons, and use the function Ws(u)/u which is given as a tabulation of coefficients fitting simple algebraic terms to solution of the complete Fokker-Planck equation in two velocity dimensions (v^ and v||).

The key quantity is Ws(u), the energy (normalized to me vr2/2) imparted to the electric field Edc by an electron as it slows down. The local electric field Edc 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, Ws(u) can be represented in a power series. From the first two terms, one finds

 Jrf = -e neG óõ dv||Dql(v||) ¶fe(v||)¶v|| v||3  45 + Z éê ë m- 1 + Z/2 + 3m2/23+Z v||2vr2 ùú û ,
(34)
where m = -1 for cooperative current drive, and m = +1 for current drive into an opposing field.

The function Ws(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 Edc 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 rf-driven 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 (n-1) the local total current density is J(n-1), the rf-driven current density is Jrf(n-1), and the local electric field is Edc(n-1). We want to advance quantities to the new time, (n), and be consistent with Ohm's Law: J = sEdc + Jrf, where s is the parallel neoclassical conductivity. An initial guess at what the new electric field Edc(n-) must be is formed from the new total current, as advanced by TSC, and the old rf-driven current and old electric field.

 J(n) = sEdc(n-)+Jrf(n-1)(Edc(n-1))     ,
(35)
From the Edc(n-) obtained from equation (35) one can find the estimated new rf current density Jrf(n-) from equation (33). Its derivative with respect to Edc,  ( Jrf/ Edc)(n-), is formed by numerical differentiation. Algebra yields estimates for the new electric field and for the new rf-driven current density,

 Edc(n) = Edc(n-) + Jrf(n-1) - Jrf(n-) s+ ( ¶Jrf/ ¶Edc)(n-) ,
(36)

 Jrf(n) = Jrf(n-) + ( ¶Jrf/ ¶Edc )(n-)( Jrf(n-1) - Jrf(n-) ) s+ ( ¶Jrf/ ¶Edc)(n-) .
(37)

The new values Edc(n) and Jrf(n) may not be consistent. This is checked by assigning them to the previous step labeled (n-1) 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

 tr-loss .n er

ne
= tr-lossDql(vr) ê
ê
ê
fe(vr)
v||
ê
ê
ê
<<  1
(38)
where tr-loss is a loss time for electrons pushed by rf waves into a runaway region and [n\dot]er is rate of increase of the density of runaway electrons owing to the rf diffusion from a lower velocity. A simple estimate for tr-loss would be the confinement time ( ~ 10 msec). If equation (38) is not satisfied, then the calculation is inappropriate.

Iteration is not necessary if when the ambient electric field is low and the density is high, such that w << |vr k|||. In that case the evolution of the electric field and the current is handled without special attention as the transport calculation evolves.

### 2.5  Broadening of rf-Driven Current

The rf-driven current is subject to diffusion across the confining magnetic field. Magnetic turbulence is often proposed as the dominant mechanism.

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, RF-driven current Jd is modified from the computed RF-driven current Jrf by slowing-down, characterized by an inverse-time nslow, and a cross-field diffusion, characterized by DJrf.

Assuming for the moment slab geometry with x the spatial variable the diffusion-like equation we use is:

 ¶Jd/ ¶t = nslow ( Jrf- Jd)  +  ¶/¶x  ( DJrf ¶Jd/ ¶x )
(39)

This can be transformed to a normalized flux coordinate [^x] º (y- ymin)/(ymax - ymin) with standard transformations for curvilinear coordinates. Under the reasonably good assumption that y ~ r2 the result is:

 ¶Jd/ ¶t = nslow ( Jrf- Jd)  +  (4/a2)  ¶/¶ ^x ( ^x DJrf ¶Jd/ ¶ ^x )    ,
(40)
where a represents the nominal minor radius of the plasma.

Equation 40 clearly has desirable properties: if DJrf is small then the Jd moves to Jrf on the nslow time scale; and if DJrf is large then Jd and Jrf can be quite different, in particular Jd will be smooth compared to Jrf over a length of order Ö{ DJrf /  nslow}.

In LSC DJrf on input is one number, constant across the plasma, in units of meter2 per second. Internally, DJrf is an array, so generalization would be simple. The nslow 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 Jd comes from one inversion of a tridiagonal matrix.

The boundary condition on Jd 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 dJd/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 Jd and dJd/d[^x] at the first interior grid point.

The rf-driven current found Jrf leads directly to the smoothed (diffused) Jd after the inversion of a tri-diagonal matrix [41].

Note that the integral of Jd over the plasma cross section, or total diffused rf-driven current, is not constrained to have any particular relationship with the same integral over Jrf, 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 mid-radius and mildly peaked density profile.

The heuristic smoothing method of V. Fuchs [17] is similar to the method described here, but apparently treats DJrf and nslow as constants, sets the change in total current from an estimate incorporating the formula for nslow, and finds the value needed for DJrf 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.

### 2.6  Heuristic Power Diffusion Estimate

The undiffused current Jrf was computed from the ray tracing, quasi-linear development of an electron distribution function in the parallel velocity giving the local power deposition per unit volume P0, and finally the Karney-Fisch [12] electric field adjustment to current density.

As explained in the previous section, the diffused current density Jd can be quite different from Jrf, and therefore quite different in radial distribution from the P0 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

P( ^
x

) = as
 |Jd( ^x )|  ne( ^x )

 óõ |Jd( ^x )|  ne( ^x )  dV( ^x )
ó
õ
P0 dV  + (1-as) P0( ^
x

)     ,
(41)
where P is the diffused power, P0 is the power deposited from the ray tracing and quasilinear calculation, as ranges from 0 (no spreading of P0) to 1 (full spreading), dV is a volume element.

### 2.7  X-ray Camera Image

The PBX-M experiment to which we have applied our calculation the most has a 2-dimensional x-ray camera with an axis roughly tangent to a mid-plane 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 2-d image with the following formula.

 Ix-p
 =
 1dAp óõ dl dA(l) dWp(l)  ni(y(l))  ne(y(l)) ×
 óõ dv||v||fe(v||,y( l)) ×
 óõ E(v||) 0 dE¢E¢ ¶2 sei¶E¢¶W [E(v||), E¢, ^n · ^b ] T(E¢)
(42)

In the above, Ix-p is the x-ray intensity on pixel p of area dAp 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, dWp is the solid angle subtended by the pinhole as seen from the emitting volume element, ni(e) is the ion (electron) density at the emitting location, E¢ denotes photon energy emitted by an electron of energy E, 2 sei/E¢W is the cross section for electron-ion 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 fe is that computed from the one-dimensional Fokker-Planck model, and is itself not influenced by the electric field Edc. The absorption from vessel and windows is computed from fits to experimental absorption data [43].

## 3  LSC Usage

The LSC code is written in standard Fortran and has run on all Unix/Linux machines available at PPPL, including the Cray machines at NERSC , and VAXes. There are no mathematical library calls; interpolators, integrators, root solvers, look-up searchers, and Bessel function approximations are part of the source code. Routines are taken directly from Numerical Recipes [41] are: LAGUERnr, ZROOTSnr, RAN3nr, HUNTnr, PIKSRTnr, and PIKSR2nr, and TRIDAInr. The trailing nr'' denotes the source, where as the preceding characters are the same as in the book [41].

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, MS-DOS, 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 3-D-like 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 2-D 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 no-operations 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 stand-alone 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 stand-alone use with a simple driver program is easily done. The distribution provides a driver, lscdrive.F , and a simplifed mock-equilibrium writer, WrCircEq.F . Stand-alone runs can be extremely useful in studying the behavior and assuring reasonable convergence with the parameters chosen.

Naturally, stand-alone computations cannot show the interaction of the current driven with the electric field generated, reacting to the current, in a complete calculation.

### 3.1  Input Files

Call to LSC, whether from a transport code or a stand-alone driver, needs to specify

1. The plasma (geometry, fields, density, temperature); This specification is done through a large, strictly formatted, equilibrium file. The name given this file is arbitrary, but stand-alone runs have commonly used jardin.d and TSC runs have commonly used lhcdoua . For TRANSP the file'' is in memory; in other words there are no disk read/writes for the equilibrium. (Im vague on this; needs checking.)
2. LH parameters such as the frequency, the spectrum, and current/power diffusion. This specification is done through short NAMELIST (list-directed) file input.lhh under the namelist name inpvalue
3. Computational paramers such as the step size and the number of steps for each ray, strategy for smoothing in velocity space , number of iterations in constructing the distribution function. This specification is done through the same short NAMELIST file input.lhh under the namelist name inpexprt . (Putting the two namelists in different files might be a good idea in some cases, and doing so would be simple enough.)
4. Flags to specify graphs and printouts. This specification is also under the namelist name inpexprt .
5. LHCD power, i/o unit numbers, a flag on the strategy for re-tracing rays, and a flag which can turn off all plots . This specification comes in the argument list in the statement calling LSC.

Details follow in the following sub-sections.

#### 3.1.1  Plasma

The data on the plasma (geometry, fields, density, temperature) is read by a subroutine rdTSC in the source file grapgrd.F .

A header line (informational comments) equhd is read first followed by integers npsitm, nspc, kcycle, a plasma-time 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.

#### 3.1.2  LH parameters

The specification is done through the NAMELIST (list-directed) 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.

Table 1: LSC input parameters related to the launched lower hybrid rays

 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 Gaussian-model 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 Gaussian-model 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 Maxwellian-model 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 nq 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.

Table 2: LSC input parameters related to lower hybrid rays computed via the Brambilla method

 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 eight-character 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 ith guide in cm; bz(i) is the toroidal distance to the far edge of ith guide in cm; so that az(2)-bz(1) is the the width of the septum between the first two guides.

Table 3: LSC input parameters related diffusing the deposited power and current

 DiffuJrf The heuristic diffusivity for rf-driven 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.

#### 3.1.3  Computational parameters

Computational paramers such as the step size and the number of steps for each ray, strategy for smoothing in velocity space , number of iterations in constructing the distribution function. This specification is done through the same short NAMELIST file input.lhh under the namelist name inpexprt. (Putting the two namelists in different files might be a good idea in some cases, and doing so would be simple enough.)

Table 4 gives the computational parameters, and lists the defaults for each one.

Table 4: LSC computational input parameters

 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 Dql 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 fe, Dql and Prf; some of these steps may be with a constant power (controlled by nFlat, above). The default is NRAMPDIM (200). nFlat The number of ramp-up 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 ramp-up iteration. The default is 0.20 - 20 percent new fe and 80 percent old fe. 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 negatively-directed'' 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.

#### 3.1.4  Flags for graphs and print files

Table 5: LSC flags for plotting and printing: flag number and internal name

 Plot Flag PlFlg(number): Descripton RAYPL 1: Ray position in r, z, f and in k^,k||; enhancement SPECPL 2: Launched ray power spectrum vs n|| and v RFDPL 3: Pray, Jray at 12 y, 6 per page; total Pray Jray vs n|| and v|| RFDPSPL 4: Pray Pql Jrf and integrals vs root y DAMPL 5: quasi-linear absorption vs root y, 6 rays per page JRFPL 6: ne Te Edc Itsc dJ/dE (E/J) PITPRFPL 7: Bz/Bf, pitch, n|| enhance, ne Prf vs Rmajor DQLFEPL 8: Dql fe vs v|| at several y Write Flag PrFlg(number): Description RAYWR 1: much ray information in labeled columns FSTFRCWR 2: fractions of fast particles NPAPWRWR 3: n||'s and powers in rays vs index

Table 6: LSC input parameters controling output

 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 ramp-up iterations the results of fe and Dql 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.

Table 7: Typical LSC namelist inputs

&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,
/

#### 3.1.5  Calling arguments

Table 8 lists the arguments that must be supplied by the calling program, and is largely self-explanatory. 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 self-consistent way iRayTrsi is crucial, and almost always equal to 0 because of the number of iterations needed.

Table 8: LSC subroutine arguments, to be supplied by calling program

 LhPwrMW Lower Hybrid power in MW passed from calling program, TSC/TRANSP/etc nTSCwrit Unit number to which caller writes equilibrium data nTSCread Unit number from which caller reads Prf, Jrf , etc nTSCscrn Unit number caller uses for writing to the screen nTSCgraf Unit number to which LSC commands and data for gnuplot; the name used is date-time.gnp nLSCcomm Unit number to which LSC sends screen output nTSCunus Unit number not used by caller; used internally by LSC iRayTrsi 0 = use old ray data, old fe(v), and use new Edc for the current; 1 = calculate new rays and fe(v) from new equilibrium; 2 = use old ray data, but calculate new fe(v) taking account of new ne and Te and use new Edc for the current iPlotsig 0 = do not make plots; 1 = make plot files asked in input.lhh iError 0 = LSC finishes without errors; 1 = one or more error(s) found; -1 = LSC found an error; calling program can keep going

### 3.2  Output files

Table 9: LSC output files

 Wr2TSC.out Electron power, electron current, E/J   dJ/dE, ion power LSCcomm.out LSCcom2.out ray.dat

### 3.3  Output graphs

Figure 1: Graphed'' input parameters to LSC for keeping track of the nature of the calculation by putting important information in the graphics output file.

Figure 2: A launched spectrum from a Brambilla calculation, including the negative-going minor components. Upper-right: Bar graph showing relative powers versus v/c . Lower-left: Same, but represented versus n-parallel. Lower-right: The launched spectrum as smoothed in velocity by the same smoothing function used in constructing the quasi-linear distribution function.

Figure 3: A typical result of the calculation for one ray. Left side: Evolution of n-parallel (top) and n-perpendicular (bottom) versus square root of normalized poloidal flux. The dotted line on the top frame indicates the region of strong damping, ie, if the magnitude of n-parallel is larger than the magnitude shown dotted, then the linear damping is strong. Quasi-linear damping may not be strong, owing to quasi-linear burn-through. Right side: Ray trajectory projected onto a poloidal cross section of the plasma. The roundness'' comes from an idealized test case; in general the cross section is elongated with the correct Shafranov shift of flux surfaces.

Figure 4: Power and current deposited without diffusing the current or the power. More rays would produce a smoother result. Horizontal axis is the square root of normalized poloidal flux. Left two: Volumetric and integrated power from the ray point of view for 1 MW launched power. Center two: As on left, but from the quasi-linear power deposition point of view. Right two: Volumetric and integrated current. The negative current stems from negative-going components of a realistic launched spectrum. The power and current found represented volumetrically and integrated from the center outwards.

Figure 5: Power and current deposited with diffusing the current and then the power. The current diffusion coefficient was set to 0.05 meter-squared per second, and the weighting of the diffusion effect on the power was 90 percent. Note that the center frames show un-diffused power.

Figure 6: A typical electron distribution and quasilinear diffusion coefficient

## 4  Appendix: Brief Description of Source Files

The following is a list of the FORTRAN source files for LSC:

Table 10: Names and functions of fortran files for LSC

 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 fe; do smoothing and derivative fits.F Fits to current response from the Karney-Fisch 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 stand-alone 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 Dql

Table 11: Names and functions of include files for LSC

 Name of Module Brief Description of Role in the Computation CGSetc.inc Constants in the cgs system Doflags.inc Flags for do-options, such as DoBram DqlBins.inc Quantities for Dql Escan.inc User-input Edc for stand-alone testing FeBins.inc Quantities for fe Implic.inc IMPLICIT NONE' and #define s for CPP Inpval.inc Namelist inputs Jrf.inc Quantities for Jrf 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 ramp-up 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

Table 12: Names and functions of include files for LSC for the X ray camera

 Name of Module Brief Description of Role in the Computation Xray.inc Quantities describing x-ray 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 x-ray camera xgraphs.inc Working arrays xparams.inc Parameters for computational bins

Table 13: Input files for LSC

 Name of File Brief Description of Role in the Computation input.lhh Parameters for the computation jardin.d Equilibrium plasma information for test cases lhcdoua Equilibrium plasma information written by TSC ray.dat Information on rays from a previous LSC run input.xry Parameters for the X-Ray camera (optional)

Table 14: Output files for LSC

 Name of File Brief Description of Role in the Computation Wr2TSC.out Power deposition and RF-driven current ray.dat Information on rays from the present LSC run date-time.gnp Combined commands and data for gnuplot

Brief Description the X-ray Camera

Table 15: Typical namelist inputs for the X-ray camera

&inpxry
E_max = 0.20,  E_min = 0.01,  E_ph_min = 1.0e-6,
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,
/

Table 16: Variables used in the X-ray camera: I

 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 two-CHARACTER parameter specifying the type of metal foil that the X-Rays 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.

Table 17: Variables used in the X-ray camera: II

 R/z_tangent The pinhole aperture of the X-Ray 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 z-position (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 X-Ray 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 X-rays by the foil when calculating the X-ray camera image; the default value of iAbsXrays is 1 (absorption). nEbins This INTEGER specifies the number of X-Ray 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 X-Ray camera image, in the x- and y-directions, 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 cross-section is divided into, radially and height-wise. 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 X-Ray 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.

## Acknowledgments

This user's manual and the development of the computer program LSC was funded by the United States Department of Energy , Office of Fusion Energy , under contract DE-AC02-76-CH03073. The computer code and the interfacing with TSC was done by D. W. Ignat , E. J. Valeo , and S. C. Jardin . Interface with TRANSP was done by D. W. Ignat, T. B. Terpstra , D. C. McCune , and S. M. Kaye .

F. W. Perkins devised the approach taken in the consideration of scattering. M. Ono , F. Paoletti , and H. Takahashi contributed discussions regarding the limited up-shift of parallel index of refraction. S. von Goeler , S. E. Jones , and D. P. Enright provided valuable help regarding the x-ray 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 re-uses many of his contributions.

## References

[1]
M. Brambilla, Slow-wave launching at the lower hybrid frequency,'' Nucl. Fusion 16 47 (1976).

[2]
S. Weinberg, Eikonal method in magnetohydrodynamics,'' Phys. Rev. 126 1899 (1962).

[3]
I. B. Bernstein, Geometric optics in space- and time- varying plasmas,'' Phys. Fluids 18 320 (1975).

[4]
J. M. Wersinger, E. Ott, and J. M. Finn, Ergodic behavior of lower hybrid decay wave trajectories,'' Phys. Fluids 21 2263 (1978).

[5]
Yu. F. Baranov and V. I. Fedorov, Lower hybrid wave propagation in tokamaks,'' Soviet Physics Technical Physics Letters 4 332 (1978); Nucl. Fusion 20 1111 (1980).

[6]
P. T. Bonoli and E. Ott, Accessibility and energy deposition of lower hybrid waves in a tokamak with density fluctuations,'' Phys. Rev. Lett. 46 424 (1981).

[7]
P. T. Bonoli and E. Ott, Toroidal and scattering effects on lower hybrid propagation,'' Phys. Fluids 25 359 (1982).

[8]
D. W. Ignat, Toroidal effects on propagation, damping, and linear mode conversion of lower hybrid waves,'' Phys. Fluids 24 1110 (1981).

[9]
F. C. Jobes, et al., Current ramp-up by lower hybrid waves in the PLT tokamak,'' Phys. Rev. Lett. 55 1295 (1985).

[10]
C. F. F. Karney, N. J. Fisch, and F. C. Jobes, Comparison of the theory and practice of lower hybrid current drive,'' Phys. Rev. A 32 2554 (1985).

[11]
C. F. F. Karney and N. J. Fisch, Numerical studies of current generation by rf traveling waves,'' Phys. Fluids 22 1817 (1979).

[12]
C. F. F. Karney and N. J. Fisch, Current in wave-driven plasmas,'' Phys. Fluids 29 180 (1986).

[13]
N. J. Fisch, Theory of current drive in plasmas,'' Rev. Mod. Phys. 59 175 (1987).

[14]
P. T. Bonoli and R. C. Englade, Simulation model for lower hybrid current drive,'' Phys. Fluids 29 2937 (1986).

[15]
E. J. Valeo and D. C. Eder, Numerical modeling of lower hybrid heating and current drive,'' J. Comput. Phys. 69 341 (1987).

[16]
P. T. Bonoli, M. Porkolab, Y. Takase, and S. F. Knowlton, Numerical Modeling of the lower hybrid heating and current drive experiments in Alcator C,'' Nucl. Fusion 28 991 (1988).

[17]
V. Fuchs, I. P. Shkarofsky, R. A. Cairns, and P. T. Bonoli, Simulations of lower hybrid current drive and ohmic transformer recharge,'' Nucl. Fusion 29 1479 (1989).

[18]
K. Kupfer and D. Moreau, Wave chaos and the dependence of LHCD efficiency on temperature,'' Nucl. Fusion 32 1845 (1992).

[19]
H. Takahashi, et al., Accessibility for lower hybrid waves in PBX-M,'' in 1993 Conference on Controlled Fusion and Plasma Physics (Proc. Conf. Lisbon, 1994) Vol 17C, Part III, European Physical Society, Geneva (1993) 901.

[20]
J. Kesner, et al., Simulations of lower hybrid current drive in shaped tokamaks,'' Nucl. Fusion 34 619 (1994).

[21]
F. Paoletti, et al., LHCD Accessibility Study with Reconstructed Magnetic Equilibria in PBX-M,'' Nucl. Fusion 34 771-776 (1994).

[22]
K. Kupfer, et al., Statistical theory of wave propagation and multi-pass absorption for current drive in tokamaks,'' Phys. Fluids B 5 4391 (1993).

[23]
R. Hawryluk, in Physics of Plasmas Close to Thermonuclear Conditions, (Varenna, 1979), Commission of the European Communities, Brussels, I 61 (1979).

[24]
R. J. Goldston, D. C. McCune, H. H. Towner, et al., J. Comput. Phys. 43 61 (1981).

[25]
S. C. Jardin, N. Pomphrey, and J. DeLucia, Dynamic modeling of transport and positional control of tokamaks,'' J. Comput. Phys. 66 481 (1986).

[26]
W. M. Nevins, et al., Mission and Design of the Tokamak Physics Experiment / Steady State Advanced Tokamak,'' Fourteenth International Conference on Plasma Physics and Controlled Nuclear Fusion Research, Würzburg, Germany, 30 September - 7 October 1992, paper IAEA-CN-56/F1-5.

[27]
J. A. Schmidt, et al., The Design of the Tokamak Physics Experiment,'' J. Fusion Energy 12 455 (1993).

[28]
K. Tani, M. Asumi, and R. S. Devoto, Numerical analysis of 2D MHD equilibrium with non-inductive plasma current in tokamaks,'' J. Comput. Phys. 98 332 (1992).

[29]
R. S. Devoto, et al., Modeling of lower hybrid current drive in self-consistent elongated tokamak equilibria,'' Nucl. Fusion 32 773 (1992).

[30]
J. P. Bizarro and D. Moreau, On ray stochasticity during lower hybrid current drive in tokamaks,'' Phys. Fluids B 5 1227 (1993).

[31]
J. P. Bizarro, On the dynamics of the launched power spectrum during lower hybrid current drive in tokamaks,'' Nucl. Fusion 33 831 (1993).

[32]
D. W. Ignat, E. J. Valeo, and S. C. Jardin, Dynamic Modeling of Lower Hybrid Current Drive,'' Nucl. Fusion 34 837-852 (1994).

[33]
S. von Goeler, et al., Initial lower hybrid current drive results and first two-dimensional images of hard x-ray emission in PBX-M,'' in Proceedings of the Nineteenth European Physical Society Conference on Controlled Fusion and Plasma Physics, (Innsbruck, Austria, 1992) II 949 (1992).

[34]
S. E. Jones, et al., Calculation of an upper limit for an effective fast electron diffusion constant using the hard x-ray camera on PBX-M,'' Plasma Phys. Control. Fusion 35 1003 (1993).

[35]
J. P. Bizarro, et al., On self-consistent ray-tracing and Fokker-Planck/ modeling of the hard x-ray emission during lower-hybrid current drive in tokamaks,'' Phys. Fluids B 5 3276 (1993).

[36]
S. von Goeler, et al., A camera for imaging hard x-rays from suprathermal electrons during lower hybrid current drive on PBX-M,'', Rev. Sci. Instrum. 65 1621 (1994).

[37]
D. W. Ignat, R. Kaita, S. C. Jardin, and M. Okabayashi, Spreading of Lower Hybrid Wave Driven Currents in PBX-M,'' Nucl. Fusion 36 1733-1742 (1996)

[38]
Myunghee Ju, Jinyong Kim, KSTAR Team, A predictive study of non-inductive current driven KSTAR tokamak discharge modes using a new transport-heating simulation package,'' Nucl. Fusion 40 1859-1866 (2000)

[39]
W. Daughton, TSC Users Manual,'' internal PPPL working document dated August 1993; updated January 1997 by W. Daughton, S. C. Jardin, C. Kessel, N. Pomphrey

[40]
S. C. Jardin, TSC Users Manual,'' 1999.

[41]
William H. Press, Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling, Numerical Recipes; The Art of Scientific Computing - FORTRAN Version, (Cambridge University Press, New York, 1989).

[42]
H. W. Koch and J. W. Motz, Bremsstrahlung Cross-Section Formulas and Related Data,'' Rev. Mod. Phys. 31 920 (1959).

[43]
W. H. McMaster, et al., `Compilation of X-ray Cross Sections,'' University of California, Livermore, California, Report UCRL 50174, Sec. II rev. 1, May 1969.

[44]
T. H. Stix, Waves in Plasmas, (American Institute of Physics, New York, 1992).

## Index (showing section)

AcDc.F, 4-0
achrg, 3-1
amass, 3-1
anecc, 3-1
anicc, 3-1
anti-Hermitian, 2-2
apl, 3-1
bgzero, 3-1
brambJES.F, 4-0
camera.inc, 4-0
centers(), 3-1
CGSetc.inc, 4-0
code
location of source, 1-0, 3-0
coordinate system, 2-1
Coulomb collision, 2-3
couplers(), 3-1
current diffusion, 2-5
current driven, 2-4
cycle.F, 4-0
cyclotron frequency
electron, 2-2
ion, 2-2
DAMPL, 3-1
date-time.gnp, 3-1, 4-0
dielec.inc, 4-0
dielectric tensor elements, 2-2
DiffuJrf, 3-1
dispersion relation, 2-2
distribution function, 2-2
Do1Rpr, 3-1
DoBram, 3-1
Doflags.inc, 4-0
DqlBins.inc, 4-0
DQLFEPL, 3-1
electric field, 2-4
emitter.inc, 4-0
emparams.inc, 4-0
Escan.inc, 4-0
ezsg.F, 4-0
fe.F, 4-0
FeBins.inc, 4-0
fghz, 3-1
fits.F, 4-0
FSTFRCWR, 3-1
g77, 3-0
gary, 3-1
gnuI.inc, 4-0
gnup.F, 4-0
gnuplot, 3-0
gnuR.inc, 4-0
gpary, 3-1
grap.F, 4-0
grapgrd.F, 3-1, 4-0
gridgen.F, 4-0
grids.F, 4-0
Hermitian, 2-2
HstpLH, 3-1
i/o unit numbers, 3-1
idiag(), 3-1
iError, 3-1
Implic.inc, 4-0
inpexprt, 3-1
inpexprt, 3-1
input files
input.lhh, 3-1
jardin.d, 3-1
lhcdoua, 3-1
input.lhh, 4-0
input.xry, 4-0
jardin.d, 4-0
lhcdoua, 4-0
ray.dat, 4-0
input.lhh, 3-1, 4-0
input.lhh, 3-1
input.xry, 4-0
Inpval.inc, 3-1, 4-0
inpvalue, 3-1
iplim, 3-1
iPlotsig, 3-1
iRayTrsi, 3-1
isym, 3-1
iteration
E-field - current, 2-4
electron distribution, 2-3
quasilinear diffusion, 2-3
jardin.d, 4-0
jardin.d, 3-1
jrf.F, 4-0
Jrf.inc, 4-0
JRFPL, 3-1
kcycle, 3-1
Landau damping, 2-2
lhcdoua, 4-0
lhcdoua, 3-1
LhPwrMW, 3-1
lower hybrid approximation, 2-2
lower hybrid resonance, 2-2
lower hybrid waves, 2-2
LSC code
location of source, 1-0, 3-0
LSC manual, location of, 1-0, 3-0
LSCcom2.out, 3-2
LSCcomm.out, 3-2
lscdrive.F, 3-0, 4-0
Makefile, 3-0
manual, location of, 1-0, 3-0
matr.F, 4-0
MKSetc.inc, 4-0
module, use as, 3-0
namelist input
inpvalue, 3-1
inpexprt, 3-1
plotting flags, 3-1
NERSC, 3-0
nFlat, 3-1
nfreq, 3-1
nGrps, 3-1
nLSCcomm, 3-1
NPAPWRWR, 3-1
nparmax, nparmin, 3-1
npolmax, npolmin, 3-1
npols, 3-1
npsi, 3-1
npsitm, 3-1
nRampUp, 3-1
nrays, 3-1
nslices, 3-1
nsmoo, 3-1
nsmw, 3-1
nspc, 3-1
nstep, 3-1
ntors, 3-1
nTSCgraf, 3-1
nTSCscrn, 3-1
nTSCunus, 3-1
nTSCwrit, 3-1
Numerical Recipes, 3-0
HUNT, 3-0
LAGUER, 3-0
PIKSRT, 3-0
PIKSRT2, 3-0
RAN3, 3-0
TRIDIA, 3-0
ZROOTS, 3-0
numerics.inc, 4-0
nv, 3-1
nx, 3-1
nz, 3-1
nzones, 3-1
output files, 4-0
date-time.gnp, 4-0
ray.dat, 4-0
Wr2TSC.out, 4-0
params.inc, 3-1, 4-0
pary, 3-1
pbxio.F, 4-0
permittivity of free space, 2-2
phaseDeg(), 3-1
PIetc.inc, 4-0
PITPRFPL, 3-1
plasejv.F, 4-0
plasma frequency
electron, 2-2
ion, 2-2
plcmx.inc, 4-0
PlFlg(), 3-1
plot flag, 3-1
plotting flags, 3-1
PlPr.inc, 4-0
poloidal flux, 2-1
power, 2-2, 2-3
power level, 3-1
power.F, 4-0
power.inc, 4-0
powers(), 3-1
ppary, 3-1
PrFlg(), 3-1
PrfSpred, 3-1
ProfBody.inc, 4-0
psilim, 3-1
psimin, 3-1
ql.F, 4-0
quasilinear diffusion, 2-3
Ramppwr.inc, 4-0
ray.dat, 4-0
ray.dat, 3-2
RayBins.inc, 4-0
Rayini.F, 4-0
Rayio.F, 4-0
RAYPL, 3-1
Raystore.F, 4-0
Raytrace.F, 4-0
RAYWR, 3-1
RayWrk.inc, 4-0
rdTSC, 3-1
re-tracing strategy, 3-1
RFDPL, 3-1
RFDPSPL, 3-1
rgzero, 3-1
rho, 3-1
runaway probability, 2-4
ScatKdeg, 3-1
sgIbk.inc, 4-0
SGlib (Scientific Graphics), 3-0
sgRbk.inc, 4-0
shells, 2-2
smoothing, 2-3
SPECPL, 3-1
spectrum, 2-2
stand-alone, use as, 3-0
Stix, T. H., 2-2, 2-3
tekev, 3-1
thermal velocity
electron, 2-2
ion, 2-2
tikev, 3-1
Tokamak Simulation Code (TSC), 1-0
TRANSP, 1-0
TSC (Tokamak Simulation Code), 1-0
TSC Users Manual, 1-0
TSCgrap.inc, 4-0
tscunits.inc, 4-0
TurnNegs, 3-1
under-relaxation, 2-3
voltlp, 3-1
vptemp, 3-1
wave diffusion, 2-3
WeghtItr, 3-1
widths(), 3-1
WkAry.inc, 4-0
WKB approximation, 2-2
Wr2TSC.F, 4-0
Wr2TSC.out, 4-0
Wr2TSC.out, 3-2
WrCircEq.F, 3-0, 4-0
write flag, 3-1
x-ray camera, 2-7
xgraphs.inc, 4-0
xmag, 3-1
xparams.inc, 4-0
Xray.inc, 4-0
XrayCam2.F, 4-0
xsv, 3-1
zmag, 3-1

File translated from TEX by TTH, version 2.79.
On 19 Nov 2000, 14:47.