DMC 14 Aug 2000
updated 24 Oct 2008
updated 05 Jan 2009, M. Gorelenkova(MG)

*******MG
The distribution function calculation is independent of the NLFBMFLR switch,
and have the code accumulate two copies of the distribution function, 
both the particle version (FBM_PTCL) and the guiding center version (FBM) .
We will still use NLFBMFLR, but only to control which copy of the 
distribution function is used to compute the beam-beam neutrons 
in ../nubeam/bbneut.for and associated routines.
*******MG


The fortran program `get_fbm' has been ported to unix and upgraded
to allow greater flexibility of access to fast ion distributions.  
In particular, beam and fusion product distribution functions (now),
and RF tail ion distribution functions (eventually we hope) can be 
accessed, at any radial or poloidal location.

Averaging and integration features allow reconstruction of such 
quantities as the flux surface averaged total fast ion density 
as a function of flux surface label, or, at a particular flux 
surface, the poloidal variation of the fast ion density.  Or the 
density integration can be restricted to particles within a certain
energy or vpll/v range.  Etc., etc.

The `get_fbm' program can also write the distribution function data
out to a NetCDF file, which will be a more convenient format for 
separately written software using this data, than the "encoded ascii"
legacy format used in the .DATAn files output directly by TRANSP.

** a note on units & normalization **

`get_fbm' reads the TRANSP COMMON array fbm from an ACfile (see
Heidbrink's note, below).  The units of FBM are

 #/cm3/eV/[delta(solid-angle)/4pi]

but in TRANSP & NUBEAM, "delta(solid-angle)/4pi" is equivalent to

   1/[number of vpll/v bins].

------------------------------------------------------------------------
(Discussion of delta(solid-angle) on a unit sphere in velocity space)

In terms of angles on a unit sphere in velocity space:

   vpll/v = cos(alpha)
   vperp/v = sin(alpha)
   delta(vpll/v) = delta(cos(alpha))

For contribution to surface area on the sphare, there is a factor 
2*pi*sin(alpha) = 2*pi*(vperp/v) since vperp/v is the radius of a 
circle on the unit sphere in 3d v-space.  But the differential area 
on the surface of the sphere also has a factor of differential width 
dw: dA=dw*2*pi*sin(alpha).  Corresponding to delta(cos(alpha)) this 
width dw of the differential band around the circle at radius 
2*pi*sin(alpha) has a factor of 1/(sin(alpha)) which cancels out the
factor due to the radius of the circle.  You can convince yourself
of this with a picture of the sphere showing the relation of dw to
delta(cos(alpha))=delta(vpll/v).

You get integral[(vpll/v = -1 to 1){d(vpll/v)*2pi}] = 4pi = surface 
area of unit sphere.  The contribution of each pitch zone delta(vpll/v)
is equal at delta(solid angle)/4pi = (1/nznbma) in terms of 
TRANSP/NUBEAM's FBM data.

(End of discussion of delta(solid-angle)).
------------------------------------------------------------------------

Given the reconstructed TRANSP COMMON arrays (e.g. in get_fbm.for), the
number of particles of fast specie 'js' in spatial zone 'jz', is given
by:

!
! nznbma = # of vpll/v bins
! nznbmea(js) = # of energy bins 
!   (in general, different fast species have different energy grids)
! bmvol(iz) = volume of iz'th spatial zone
!
      znsum=0.0      ! sum to contain total number of particles this zone
!
      do ia=1,nznbma
         do ie=1,nznbmea(js)
            dE=efbmb(ie+1,js)-efbmb(ie,js)
            znsum=znsum+fbm(ie,ia,iz,js)*dE*bmvol(iz)/nznbma
         enddo
      enddo

Questions:  send email to dmccune (dmccune@pppl.gov).
`get_fbm' enhancements can be requested.

--------------------------------------------------------
Note by W. Heidbrink (1997)

Date: Tue, 12 Aug 1997 10:27:29 -0700 (PDT)
From: <HEIDBRINK@GAV.GAT.COM>
Subject: Beam distribution function in TRANSP
To: GREENFIEL@GAV.GAT.COM, HEIDBRINK@GAV.GAT.COM, LAZARUS@GAV.GAT.COM,
MURAKAMI@GAV.GAT.COM, TERPSTRA@GAV.GAT.COM, PETTY@GAV.GAT.COM,
FOREST@GAV.GAT.COM, STALLARD@GAV.GAT.COM, BaylorLR@ORNL.GOV

                TRANSP FAST-ION DISTRIBUTION FUNCTION
                       W.W. Heidbrink

TRANSP calculates many quantities derived from the fast-ion distribution
function "f" but does not ordinarily save "f" itself.  Using tools
originally written by McCune and at JET, Ted Terpstra and I have created
programs to view "f" at DIII-D.

-----------------------------------------------------------------------------


1) Add lines to namelist.  You must instruct TRANSP to dump its memory
at certain selected times.  This creates a so-called ACFILE in the
appropriate result area.  Here are the lines I inserted in
71524A01TR.DAT:


!
! Dump common block:
!
!SELOUT='XXX'      !list of TRANSP common symbols to output (default=all)
SELAVG='FBM BMVOL BDENS2 EBA2PL EBA2PP' !quantities to average
AVGTIM=0.030   !averaging time (to reduce MC noise)
MTHDAVG=2      !2 for averaging after each MC timestep
OUTTIM=1.875,2.150  !up to 5 times that AC files are written


------------------------------------------------------------------------

2) Run TRANSP.  If all goes well, you'll end up with files like
$TRANSPROOT/result/D3D.91/71524A01.DATA1

$TRANSPROOT/result/D3D.91/71524A01.DATA2 
with one file for each OUTTIM you selected.

-----------------------------------------------------------------------

3) Use your chosen method for accessing the fbm data; use get_fbm to 
convert the data to NetCDF format if needed.

----------------------------------------------------------------------------

5) Physics note:  the "fbm" in TRANSP actually has a complicated
spatial index label that moves around radially and poloidally.  There
are codes that only bother to return "f" at the outer midplane.
This nearly always is sufficient to reconstruct the full distribution
function, because almost every orbit crosses the outer midplane.  
But not quite.  There are some "inner" orbits (close to the x-point
"pinch" orbit in phase space) that always have R<R_0.  These are
extremely rare for DIII-D beam ions, but could be important if we ever
pull out an enormous rf tail.

============================================================================
D. McCune:  Caution: NLFBMFLR = .TRUE. is the namelist default!!
***** See comments by MG
In NUBEAM,

NLFBMFLR = .TRUE. means: accumulate fbm(rho,theta,E,vpll/v) at particle
positions; 

NLFBMFLR = .FALSE. means: accumulate fbm(rho,theta,E,vpll/v) at the orbit
guiding centers.
*****

The default is .TRUE. because of the history of use of fbm for synthetic
diagnostics code, but it is not necessarily the best choice.  In particular,
fbm has no gyro-phase dimension, and so, the accumulation of data at 
particle locations that correspond to loss orbits for some gyro-phases (or
equivalently, some guiding center locations that could contribute) but
not for others, will be accumulated mostly from non-loss orbits, but, the
resultant anisotropy in the direction of contributing ions at (E,vpll/v) at 
the spatial location cannot be reconstructed!

For NLFBMFLR=.FALSE., an FLR step needs to be taken to get to particle
positions e.g. for charge exchange eflux diagnostics simulation based on fbm.
Existing codes may not do this.

Of course, if the finite larmor radius is small, the switch will have little 
effect.  But if it is not, e.g. in STs or other low field experiments, it is 
important to consider the proper setting of NLFBMFLR, which depends on the 
use to be made of the "fbm" data.

============================================================================
D. McCune:  Description of iregular 2d spatial grid over which the 
"fbm" data is written...

The MCgrid, or "Monte Carlo grid", is an irregular 2d spatial grid that
was developed originally for capture of 2d binned data in a Monte Carlo
calculation (NUBEAM).  It is used for fast ion distribution functions and
for certain spatially 2d profiles output by NUBEAM-- such as, beam halo
thermal neutral sources, beam-target and beam-beam fusion rates.

The grid is aligned with flux coordinates.  It is constructed of "zone
rows" that are equally spaced in rho=sqrt(Psi_tor/Psi_tor_at_bdy); each
zone row is subdivided into a different number of zones equispaced in
the equilibrium poloidal angle coordinate theta, with fewer subdivisions
for zone rows near the axis, and more for zone rows in the edge.  The
result is a set of zones all with roughly equal volume and cross sectional
area, as is desirable for consistency of Monte Carlo summation statistics.

The number of poloidal zones per zone row is linear in the zone row index:
typically, 4 in the 1st zone row adjacent to the magnetic axis, 8 in the 
next row out, 12 in the next row, etc.

There are two variants, according as the underlying MHD equilibrium is
updown symmetric, or updown asymmetric.  

As internally stored, the first row layout is this for the updown SYMMETRIC
variant:
   _____
  /  |  \
 / 2 | 1 \
 |___|___|

with theta=0 on the large major radius side; increasing to theta=pi on
the large major radius side covering the upper half of the plasma cross
section above the midplane.  The next row out would have 4 zones, the
next 6 zones, etc.

And this for the updown ASYMMETRIC variant:

   _____
  /  |  \
 / 4 | 3 \
 |___|___|
 |   |   |
 \ 1 | 2 /
  \__|__/

with theta=-pi on the lower branch on the small major radius side, theta=0
on the large major radius side, and theta=pi on the upper branch on the small
major radius side.  The next zone row contains 8 zones, the next 12 zones,
and so on.

Poloidal zones are stored contiguously, with the poloidal zone index
increasing with increasing theta coordinate, theta being oriented
counter-clockwise in the plasma cross section drawn to the right of
the machine axis of symmetry.

--------------------------------------------------------
Spatial resolution of fbm output from TRANSP

There are two relevant TRANSP namelist controls:

  NZONE_NB  (default =20) -- base radial resolution for NUBEAM
                             rule: divisible by 10
  NZONE_FB  (default =10) -- radial resolution for fbm
                             rule: NZONE_FB an even divisor of NZONE_NB

Thus for NZONE_NB=20, legal values of NZONE_FB are 5, 10, and 20.
For NZONE_NB=40, legal values of NZONE_FB are 5, 10, and 20 and 40.

The number of poloidal zones per radial row depends on the radial row index:

NZONE_FB = 10

row index           updown symmetric        updown asymmetric

  1                      2                       4
  2                      4                       8
  3                      6                      12
  4                      8                      16
  5                     10                      20
  6                     12                      24
  7                     14                      28
  8                     16                      32
  9                     18                      36
 10                     20                      40
                     -----                   -----
total                  110                     220


NZONE_FB = 20

arithmetic progression extended through 20 zone row indices.  The total
number of zones is NZONE_FB*(NZONE_FB+1) for the updown symmetric case, and
twice that for the updown asymmetric case: 420 and 840, respectively, for
NZONE_FB = 20.

Some variants of the MC grid include extra zone rows that extend beyond 
the plasma boundary.