% $Id: adj.web,v 4.9 90/06/16 07:32:04 karney Exp $

% $Log:	adj.web,v $
% Revision 4.9  90/06/16  07:32:04  karney
% More documentation.
% 
% Revision 4.8  90/06/15  13:59:56  karney
% Inserted more user documentation on how to compile and run the code.
% 
% Revision 4.7  90/05/11  17:37:26  karney
% Cleaned up use of macros to provied generic UNIX version.  This should
% probably be restructure to provide more orthogonality.
% 
% Revision 4.6  90/05/07  17:42:25  karney
% Fix typo in readline addition.
% 
% Revision 4.5  90/05/07  12:24:15  karney
% Use readline library to get terminal input on Sun
% 
% Revision 4.4  90/04/26  16:22:13  karney
% Second attempt to get version number
% 
% Revision 4.3  90/04/26  16:19:54  karney
% Include version number in title
% 
% Revision 4.2  90/04/19  14:46:35  karney
% Minor cosmetic changes following John's suggestions
% 
% Revision 4.1  90/04/19  12:11:43  karney
% Initial entry of Aug 22, 1989 version of code

\def\title{Adjoint function for a toroidal plasma; $ $Revision: 4.9 $ $}

% Change some of FWEB's formatting macros.
\def\\#1{\leavevmode\hbox{\def\${{\sl \char`\$}}\def\%{\char`\%}%
                          \it#1\/\kern.05em}} % Italic type for identifiers
\def\@@#1{\leavevmode\hbox{\def\${\char`\$}\def\%{\char`\%}%
                           \rm#1\kern.05em}} % Roman intrinsic functions
\let\Csc=\rm       % Font for comments

% Some general math macros
\def\abs#1{\left\vertbar#1\right\vertbar}    % Absolute values
\def\frac#1#2{{#1\over#2}}
\let\d=\partial
\def\pd#1#2{\mathchoice{\d#1\over\d#2}{\d#1/\d#2}{\d#1/\d#2}{\d#1/\d#2}}

% Some shortcuts for this program
\def\u{{\bf u}}\def\S{{\bf S}}
\def\subb#1#2#3{_{#1[#2]#3}}
\def\ja#1{\tilde\jmath_{#1}}
\def\jn#1#2#3{\bar\jmath\subb{#1}{#2}{#3}}
\def\qa#1{\tilde q_{#1}}
\def\s{*}
\def\lave#1{\overline{#1}}

@n

@* Introduction.  This program is supposed to calculate the Green's
function for rf-driven current in an axisymmetric tokamak.  This is the
solution to the Spitzer-H\"arm problem with the loop voltage $\Phi$ set to
$TL/Q$.  The rf-driven current density is given by
$$J_{\parallel 0}= \int d^3\u_0\, \lambda \S_0\cdot\nabla_0\chi,$$
where $\S_0$ is the bounce-averaged rf-induced flux (at the midplane).

In this program, velocities are normalized to an arbitrary reference
velocity $u_n$ (which is usually either $\sqrt{T/m}$ or $c$), times to
$\tau_n=(\Gamma^{e/e}/u_n^3)^{-1}$, current densities to $qnu_n$, fluxes to
$n/(\tau_n u_n^2)$, and $\chi$ to $qu_n\tau_n$.

(1) May 3, 1989:
Initial version with only the first Legendre harmonic running for RF
conference, May 1--3, 1989.

(2) May 29, 1989:
Include an arbitrary number Legendre harmonics.

(3) May 31, 1989:
Convert to using dynamic memory allocation.

(4) August 22, 1989:
Allow choice on boundary conditions at $u_{\rm max}$.

@ This program is written to run on VAX/VMS systems, Sun Sparcstations
(using double precision) and on Crays with cft77 (using single precision).
This is controlled by four macros |VAX|, |SUN|, |CRAY|, |UNIX|.  It should be
relatively straightforward to modify to program to run in other
environments.  The major system dependent aspect of this code is the use of
dynamic memory.  However this is easily turned off using the macros
|STATIC|.  Then all memory is taken from a single statically allocated
array of size |MEMSIZE|.

For the VAX, the prgram should be created with:

{\parskip=0pt\parindent=10em\tt
ftangle -mVAX adj\par
fortran/g\_float adj\par
link adj\par}

For the Cray, use:

{\parskip=0pt\parindent=10em\tt
ftangle -mCRAY -d  adj\par
{\rm (send {\tt adj.f} to Cray)}\par
cf77 -o adj adj.f\par}

(The {\tt -d} converts |do|/|end do| constructions into labeled |do|
loops.)

For the Sun, use:

{\parskip=0pt\parindent=10em\tt
ftangle -mSUN -=c=\#.web.c adj\par
f77 -I/usr/gnu/include adj.f adj.web.c -o adj $\backslash$\par
\hskip3em -L/usr/gnu/lib -lreadline -ltermcap\par}

If you need to use static memory, then specify {\tt -mSTATIC -m"MEMSIZE
100000"} to ftangle; this uses static memory of size 100000.

@ Input parameters.  When you run this program it asks you to input the
following parameters: |t|, |c2|, |ze|, |zi|, |eps|, |nmax|, |umax|, |imax|,
|kmax|, |dt|, |tmax|, |alpha|, |rho|, |lmax|.

|t| is the temperature in units of $m u_n^2$.  Thus if you want velocities
normalized to $\sqrt{T/m}$, set $t = 1$.  Alternatively if normalization is
to speed of light, set $t = T/(mc^2)$.

|c2| is the square of the speed of light (normalized to $u_n^2$).  If you
want to run nonrelativistic cases, then set $t=1$ and ${\it c2} = 10^{20}$.

|ze| is the electron charge state.  This is usually 1, but you can get the
Lorentz limit $Z_i\to\infty$ with $|ze|=0$ and $|zi|=1$.

|zi| is the ion charge state.

|eps| is the inverse aspect ratio $r/R$.  This is only used in the magnetic
field routine |deltab|.  This routine is currently set up to define a
circular cross-section tokamak.  $\epsilon=0$ corresponds to a uniform
magnetic field (no trapped particles).

|nmax| is the number of grid points in the $u_0$ direction.

|umax| is the maximum value of $\abs{{\bf u}_0}$.  The $u_0$ grid is evenly
spaced between $0$ and $u_{\rm max}$.

|imax| is the number of grid points in the $\theta_0$ direction.  The
$\theta_0$ grid is evenly space between $0$ and $\theta_{\rm tr}$, the
trapping angle.  $\chi({\bf u}_0)$ is odd about $\theta = \pi/2$ and is
zero for $\theta_{\rm tr} < \theta < \pi - \theta_{\rm tr}$ in the trapped
zone.

|kmax| is the number of integration points for field line integrals.  These
are only done once at the beginning of a run, so |kmax| should be set quite
large (e.g., 1000).

|dt| is the time step.  This needs to be chosen large enough to get to the
steady state quickly but not so large that the iteration becomes unstable.
A typical value is 1.

|tmax| is the number of time steps, typically 100--300.

|alpha| is a weighting factor in the time stepping algorithm.  $\alpha=0$
gives an explicit algorithm.  $\alpha=0.5$ gives Crank-Nicholson iteration.
$\alpha=0.55$ (i.e., a small amount of forward weighting) seems to be best.

|rho| governs the treatment of the boundary at $u=u_{\rm max}$.  For
$\rho=1$, we have $\chi''(u_{\rm max})=0$, and for $\rho=0$, we have
$\chi'(u_{\rm max})=0$.  The boundary has the least effect on the solution
in the interior for $\rho=1$.

|lmax| is the maximum number of Legendre harmonics.  $|lmax|=0$ corresponds
to collision off a Maxwellian.  For $\epsilon=0$, set $|lmax|=1$ (since
$\chi$ consists of the first Legendre harmonic only).  Otherwise try
$|lmax|=21$.

The program cycles through repeatedly asking for input until the end of the
input file (for terminal input this is Control-Z on the VAX or Control-D on
the Sun) or else nonpositive values are read in for |t| or |c2|.

@ Output.  This program outputs a small amount of data to the terminal plus
a file of unformatted data.

The terminal output consists of an echo of the input parameters.  Every 10
time steps, a line is written out consisting of |time|, |cond|, |abserr|,
|relerr|.  |time| is the time step.  |cond| is the present estimate of the
electrical conductivity (which this program calculates because $\chi$ is
just the Spitzer-H\"arm function).  |abserr| and |relerr| are estimates of
the absolute and relative errors in $\chi$ based on the residue.

The data written to the file consists of the input variables, the $u$ and
$\theta_0$ arrays and a two-dimensional array for $\chi$.  See the section
``Write out results to disk.''  Note that $\chi$ and $f_m$ are defined at
the mid-points of the cells i.e., at $u=|ua(n)|$, $\theta_0=|th0(i)|$ for
$0\le n<n_{\rm max}$, $0\le i<i_{\rm max}$.

@ Files.  The program talks to three different files, which we call
|input|, |output|, |adjout|.  |input| and |output| are normally terminal
input and output.  |adjout| is the output file of unformatted data (see
above).

On the VAX, the system supplies a {\tt .dat} file type to the |adjout| file
(i.e., it is normally the {\tt adjout.dat} in the current directory).  You
can redirect |input| and |output| from/to a file and rename the file
|adjout| with VAX logical names.  For example:


{\parskip=0pt\parindent=10em\tt
define/user sys\$output adj2.in\par
define/user sys\$input adj2.out\par
define/user adjout adj2.dat\par
r adj\par}

On the Sun, you can use normal file redirection for |input| and |output|.  You have to rename |adjout| explicity.  For example:

{\parskip=0pt\parindent=10em\tt
adj <adj2.in >adj2.out\par
mv adjout adj2.dat\par}

On the Cray, your can use command line arguments to specify file names.
For example:

{\parskip=0pt\parindent=10em\tt
adj input=adj2.in output=adj2.out adjout=adj2.dat\par}

@ References:

[1] B. J. Braams and C. F. F. Karney, {\it Conductivity of a relativistic
plasma}, Phys.~Fluids {\bf 1B}, 1355--1368 (1989).

[2] C. F. F. Karney, N. J. Fisch, and A. H. Reiman,
{\it Green's function for rf-driven current in a toroidal plasma},
in {\it Radio-Frequency Power in Plasmas},
Proc.\ Eighth Topical Conf. on Radio Frequency Plasma Heating,
University of California, Irvine, CA,
May 1--3, 1989,
edited by R. McWilliams,
AIP Conference Proceedings, No. 190
(AIP, New York, 1989), pp.~430--433.

[3] C. F. F. Karney and N. J. Fisch, {\it Efficiency of Current Drive by
Fast Waves}, Phys.~Fluids {\bf 28}, 116--126 (1985).

[4] C. F. F. Karney, {\it Fokker--Planck and Quasilinear Codes},
Computer Physics Reports {\bf 4}, 183--244 (1986).

@ Macros for machine dependencies.

@#if defined(CRAY)
@m CRAY 1
@m SUN 0
@m VAX 0
@m UNIX 0
@#elif defined(SUN)
@m CRAY 0
@m SUN 1
@m VAX 0
@m UNIX 0
@#elif defined(VAX)
@m CRAY 0
@m SUN 0
@m VAX 1
@m UNIX 0
@#else
@m CRAY 0
@m SUN 0
@m VAX 0
@m UNIX 1
@#endif

@ Define precision.

@#if VAX || SUN || UNIX  /* Use double precision on the VAX and Sun*/
@m real     double precision /* Default precision */
@m areal    DBLE             /* Coerce to default precision */
@m const(x) x##d0            /* Constant in default precision */
@m epsilon  1.0d-16          /* Approximate round-off error */
@m single   REAL             /* Coerce to single precision */
@m memunit  8                /* Size of a real in bytes (the VAX's unit of
                                 memory) */
@#endif

@#if CRAY /* Use single precision on the Cray */
@m real     REAL
@m areal    REAL
@m const(x) x##e0
@m epsilon  1.0e-16
@m single   REAL
@m memunit  1                /* Size of a real in words (the Cray's unit of
                                memory) */
@#endif

@ I/O units.  Also set things up to use the readline library if it's
available (to permit editing and recall of input lines).

@m outunit  6
@m diskunit 20

@#if !defined(READLINE)
@#if SUN
@m READLINE 1
@#else
@m READLINE 0
@#endif
@#endif

@#if READLINE
@m LINELEN 80
@m inunit   line(1:min(linelen,LINELEN))
@#else
@m inunit   5
@#endif

@ Label definitions.

@m Start    #:0
@m Exit     #:0

@ The Sun uses |implicit undefined(a-z)| instead of |implicit none|

@f implicit_none implicit
@#if SUN
@m implicit_none implicit undefined(a-z)
@#else
@m implicit_none implicit none
@#endif



@ Abbreviated names for some identifiers.  CFT cannot handle identifiers
longer than 8 characters.

@#if CRAY
@m maxwellian    maxwel
@m potential     pot
@m isotropic     isotrop
@m matrixinit    matinit
@m legendrej     legj
@m legendrejs    legjs
@m legendrejn    legjn
@#endif

@* The main program.


@a
      program adjoint

      implicit_none
      integer nmax,imax,kmax,tmax,lmax,tskip
      real t,c2,ze,zi,eps,umax,dt,alpha,rho
      real cpu,cpu0
      external operinit,second,runa
@#if READLINE
      character line*LINELEN
      integer readline,linelen
      external readline
@#endif

      call operinit            /* Open files */

Start: continue

      write(outunit,*)
     $     'Enter t,c2,ze,zi,eps,nmax,umax,imax,kmax,dt,tmax,alpha,rho,lmax'
@#if READLINE
      linelen=readline("> ",line)
      if (linelen<=0) then
         write(outunit,*)
         goto Exit
      end if
@#endif
      read(inunit,*,end=Exit) t,c2,ze,zi,eps,nmax,umax,imax,kmax,
     $     dt,tmax,alpha,rho,lmax
      if (t<=0 || c2<=0) goto Exit
      write(outunit,*) t,c2,ze,zi,eps,nmax,umax,imax,kmax,
     $     dt,tmax,alpha,rho,lmax
      tskip=10

      call second(cpu0)

      call runa(t,c2,ze,zi,eps,nmax,umax,imax,kmax,dt,tmax,
     $     alpha,rho,lmax,tskip)

      call second(cpu)

      write(outunit,'(1x,a,f10.2,a)') 'CPU time: ',cpu-cpu0,' secs'

      goto Start

Exit: continue

      stop
      end

      @<Functions and subroutines@>

@#if SUN || UNIX
@c
      @<C functions@>
@#endif

@ Open files for the CRAY.

@m STRING(s) #s   /* Macros get unit numbers into argument of |link|. */
@m S(s)      STRING(s)

@<Functions...@>=
@#if CRAY
      subroutine operinit

      implicit_none
      external link

      call link('input=terminal,unit'//S(inunit)//'=(input,open,text),'// @|
     $     'output=terminal,unit'//S(outunit)//'=(output,create,text),'// @|
     $     'print'//S(outunit)//','// @|
     $     'unit'//S(diskunit)//'=(adjout,create,seq)//')

      return
      end
@#endif

@ Open files for the VAX.  Define a function to get the CPU time

@<Functions...@>=
@#if VAX
      subroutine operinit

      implicit_none
      external lib$init_timer

      call lib$init_timer
      open(unit=diskunit,file='adjout',form='unformatted',status='new')

      return
      end

      subroutine second(cpu)

      implicit_none
      integer cputime
      real cpu
      external lib$stat_timer

      call lib$stat_timer(2,cputime) /* |cputime| is in 10 ms increments */
      cpu=const(0.01)*cputime /* convert to seconds */

      return
      end
@#endif

@ Open files for the SUN.  Define a function to get the CPU time

@<Functions...@>=
@#if SUN || UNIX
      subroutine operinit

      implicit_none
      integer cputime,clock
      external clock

      cputime=clock()
      open(unit=diskunit,file='adjout',form='unformatted',status='unknown')

      return
      end

      subroutine second(cpu)

      implicit_none
      integer cputime,clock
      real cpu

      cputime=clock()
      cpu=const(0.000001)*cputime /* convert to seconds */

      return
      end
@#endif

@ A readline function for the Sun.  This reads in a line with editing and
recall using the Gnu readline library.

@<C functions@>=
@#if READLINE
#include <stdio.h>
#include <strings.h>
#include <readline/readline.h>

/* Read a string with a prompt returning string in |result|.  String length
is returned as function value. */
long int readline_ (prompt, result, lprompt, lresult)
char prompt[], result[];
long int lprompt, lresult;
{
  char *cprompt, *cresult;
  cprompt = (char *)malloc (lprompt+1);
  strncpy (cprompt,prompt,lprompt);
  cprompt[lprompt] = 0;

  /* Get a line from the user. */
  cresult = readline (cprompt);
  free (cprompt);

  /* If the line has any text in it, save it on the history. */
  if (cresult && *cresult)
    add_history (cresult);

  if (cresult == (char *)NULL) {
    return (EOF);
  } else {
    strncpy (result, cresult, lresult);
    lresult = strlen (cresult);
    free (cresult);
    return (lresult);
  }
}
@#endif

@ Some C compatibility routines for the Sun.

@<C functions@>=
@#if SUN || UNIX
long clock_()			/* Interface to C clock routine */
{
  return(clock());
}
@#endif

@* Memory allocation.  All memory for arrays is allocated in |runa|.  When
tangled for static memory allocation, the memory is taken from a single
statically allocated array |memory| of size |MEMSIZE|.  For dynamic memory
allocation, memory is obtained using |lib$get_vm| or |mzalloc|.  These
return an address of a memory block.  This is associated with an array with
a pointer statement (for the Cray) or by passing the address by value down
to |runb|.  (See the definition of the macro |mem|.)

@#if !UNIX
@#if defined(STATIC)
@m STATIC 1        /* Allocate memory statically? */
@#else
@m STATIC 0
@#endif
@#else
@m STATIC 1
@#endif

@#if !defined(MEMSIZE)
@m MEMSIZE 4000
@#endif

@ Utility macros for memory allocation.  For each array used in |runb|, a
index is defined by prepending a |p| to the name of the array.  This index
is then defined using |alloc|.  Arrays which are local to a single
subroutine are allocated from a single array |work|.  |workalloc| figures
out the size of this array.

@m ptr(x)        p##x  /* A pointer into |memory| formed by prepending |p| */
@m alloc(x,y)    ptr(x)=memsize; memsize=memsize+y;
@m workalloc(y)  worksize=max(worksize,y)

@#if STATIC
@m mem(x)        memory(ptr(x))
@#else
@#if VAX
@m mem(x)        %val(baseaddr+memunit*ptr(x))
@#endif
@#if CRAY || SUN
@m mem(x)        memory(ptr(x))
@#endif
@#endif

@ Since we deal only in harmonic numbers $l$ equal to 0, 1, 3, 5, etc., we
also define macros to map between these and an index |la| equation to 0, 1,
2, 3, etc.  We also need macros to index into two non-rectangular arrays
|garr| and |harr|.

@m laf(l)   (l+1)/2      /* Map |l| onto collapsed index |la| */
@m lf(la)   max(2*la-1,0)  /* Map |la| onto expanded index |l| */

@m gind(la,la1)     ((((la+3)-1)*(la+3))/2 + (la1+1))
@m gv(la,la1)       garr(gind(la,la1))
@m hind(la,la1,la2) (((la-1)*la*(2*la-1))/6 + (la1-1)*la + la2)
@m hv(la,la1,la2)   harr(hind(la,la1,la2))

@ Allocate the necessary memory and call |runb| to do the actual
computation.

@f @<Declare memory allocation stuff@> integer
@f @<Index declarations@> integer
@f pointer integer

@<Functions...@>=
      subroutine runa(t,c2,ze,zi,eps,nmax,umax,imax,kmax,dt,tmax,
     $    alpha,rho,lmax,tskip)

      implicit_none
      integer nmax,imax,kmax,tmax,lmax,tskip
      real t,c2,ze,zi,eps,umax,dt,alpha,rho

      integer nmax1,imax1,lmaxa,lmaxa1
      @<Declare memory allocation stuff@>
      external runb

 /* |nmax1| and |imax1| are the allocated dimensions for two-dimensional
 arrays.  These might differ from |nmax| and |imax| on the Crays, where even
 strides are inefficiently handled. */
      nmax1=nmax
      imax1=imax
@#if CRAY   /* To avoid memory conflicts on the Crays */
      if (mod(nmax1,2)==1) nmax1=nmax1+1   /* Make even for accessing
               |ijy(0:nmaxa1,0:6,0:1)| by row in |potential| */
      if (mod(imax1,2)==0) imax1=imax1+1   /* Make odd for accessing
               |chi(0:imax1-1,0:nmax-1)| by row in |decomp| */
@#endif

 /* |lmaxa1| is the allocated dimension for the |la| index.  This is
 usually |lmaxa| but must be at least $1$ to ensure that dimension
 specifications such as |1:lmaxa1| are legal.  Also $H_{1,1,1}$ is needed
 even when |lmaxa=0|. */

      if (mod(lmax,2)==0) lmax=lmax-1  /* Round |lmax| down to an odd number */
      lmaxa=laf(lmax)
      lmaxa1=lmaxa
      if (lmaxa1==0) lmaxa1=1        /* But |lmaxa1| needs to be at least 1 */

      memsize=0
      worksize=0
      @<Index specifications@>  /* Define all the array pointers and
                                   figure |worksize| */
      alloc(work,worksize)

      @<Allocate memory@>

      call runb(t,c2,ze,zi,eps,nmax,umax,imax,kmax,dt,tmax,
     $     alpha,rho,lmax,tskip,
     $     nmax1,imax1,lmaxa,lmaxa1,worksize,
     $     mem(ua),mem(ug),mem(th0a),mem(th0g),mem(csth0a),mem(chi),
     $     mem(legjy),mem(chil),mem(legp),mem(duug),mem(fu),mem(dtt),
     $     mem(fm),mem(psid),mem(lam1a),mem(lam2g),mem(harr),mem(uinv),
     $     mem(th0inv),mem(resid),mem(resida),mem(work))

      @<Deallocate memory@>

      return
      end

@
@<Declare memory allocation stuff@>=
      integer memsize,worksize,ptr(work)
      @<Index declarations@>
@#if STATIC
      real memory
      dimension memory(0:MEMSIZE-1)
@#else
@#if VAX
      integer baseaddr,status,lib$get_vm,lib$free_vm
      external lib$get_vm,lib$free_vm
@#endif
@#if CRAY
      real memory
      dimension memory(0:*)
      pointer (baseaddr,memory)
      integer status,mzalloc,mzdalloc
      external mzalloc,mzdalloc
@#endif
@#if SUN
      real memory
      dimension memory(0:*)
      pointer (baseaddr,memory)
      integer status,malloc
      external malloc,free
@#endif
@#endif

@ Allocate memory required using the appropriate system calls.

@#if !STATIC
@#if VAX
@m memalloc(size,addr)   status=lib$get_vm(size,addr)
@m memdealloc(size,addr) status=lib$free_vm(size,addr)
@m bad(s)                mod(s,2)==0
@#endif
@#if CRAY
@m memalloc(size,addr)   status=mzalloc(addr,size)
@m memdealloc(size,addr) status=mzdalloc(addr,size)
@m bad(s)                s!=0
@#endif
@#if SUN
@m memalloc(size,addr)   addr=malloc(size); status=addr
@m memdealloc(size,addr) call free(size)
@m bad(s)                s==0
@#endif
@#endif

@<Allocate memory@>=
@#if STATIC
      if (memsize>MEMSIZE) then
         write(outunit,*) 'Memory needed/allocated: ', memsize, '/', MEMSIZE
         return
      end if
@#else
      memalloc(memunit*memsize,baseaddr)
      if (bad(status)) then
         write(outunit,*) 'Couldn''t allocate memory: ',memsize
         return
      end if
@#endif

@
@<Deallocate memory@>=
@#if !STATIC
      memdealloc(memunit*memsize,baseaddr)
@#endif

@* Do a case.

@f @<Global declarations@> integer

@<Functions...@>=
      subroutine runb(t,c2,ze,zi,eps,nmax,umax,imax,kmax,dt,tmax,
     $    alpha,rho,lmax,tskip,
     $    nmax1,imax1,lmaxa,lmaxa1,worksize,
     $    ua,ug,th0a,th0g,csth0a,chi,legjy,chil,legp,duug,fu,dtt,
     $    fm,psid,lam1a,lam2g,harr,uinv,th0inv,resid,resida,work)

      implicit_none

      integer worksize
      real work
      dimension work(0:worksize-1)

      @<Global declarations@>

      @<Initialization@>

      do time=0,tmax
         @<Make a step@>
      end do

      @<Write out results@>

      return
      end

@* The velocity space meshes.

@<Glob...@>=
      integer nmax,imax,nmax1,imax1,n,i
      real du,umax,ua,ug,dth0,th0max,th0a,th0g,csth0a,chi @/
 /* The velocity mesh */
      dimension ua(0:nmax-1),ug(0:nmax),th0a(0:imax-1),th0g(0:imax),
     $    csth0a(0:imax-1),chi(0:imax1-1,0:nmax-1)
      real pi

@ Define the indices for velocity space meshes in main memory of |runa|.
@<Index decl...@>=
      integer ptr(ua),ptr(ug),ptr(th0a),ptr(th0g),ptr(csth0a),ptr(chi)

@ Allocate the space needed for these arrays.
@<Index spec...@>=
      alloc(ua,nmax)
      alloc(ug,nmax+1)
      alloc(th0a,imax)
      alloc(th0g,imax+1)
      alloc(csth0a,imax)
      alloc(chi,imax1*nmax)

@ Initialize velocity meshes, $\chi$, etc.

@<Init...@>=
      du=umax/nmax
      do n=0,nmax-1
         ua(n)=(n+const(0.5))*du
         ug(n)=n*du
      end do
      ug(nmax)=umax

      pi=4*atan(const(1.0))
      th0max=atan2(const(1.0),sqrt(deltab(pi,eps))) /* $\sin^{-1}(1/\sqrt b)$*/
      dth0=th0max/imax
      do i=0,imax-1
         th0a(i)=(i+const(0.5))*dth0
         csth0a(i)=cos(th0a(i))
         th0g(i)=i*dth0
      end do
      th0g(imax)=th0max

      do n=0,nmax-1
         do i=0,imax-1
            chi(i,n)=0
         end do
      end do

@* Potentials.  The following routines calculate the potentials and their
derivatives.  First define arrays to hold Legendre functions evaluated on
the velocity mesh.

@m jl(n,a,la)  legjy(n,a,0,la)
@m yl(n,a,la)  legjy(n,a,1,la)
@m djl(n,a,la) legjy(n,a,2,la)
@m dyl(n,a,la) legjy(n,a,3,la)

@<Glob...@>=
      real c,c2
      real legjy            /* Arrays holding Legendre functions */
      integer lmax,lmaxa,lmaxa1
      dimension legjy(0:nmax1-1,0:6,0:3,0:lmaxa1)
      external initjy

@
@<Index decl...@>=
      integer ptr(legjy)

@
@<Index spec...@>=
      alloc(legjy,nmax1*7*4*(lmaxa1+1))

@ Fill the arrays.
@<Init...@>=
      c=sqrt(c2)
      lmaxa=laf(lmax)
      call initjy(ua,c,nmax,nmax1,lmaxa,lmaxa1,legjy,
     $     work(0),work(2*(lf(lmaxa)+2+1)*7),lf(lmaxa)+2)

@ Define |initjy|.
@<Functions...@>=
      subroutine initjy(ua,c,nmax,nmax1,lmaxa,lmaxa1,
     $    legjy,jarray,djarray,lmax1)

      implicit_none
      integer nmax,nmax1,lmaxa,lmaxa1,n,mult,la,l,a,lmax1
      real ua,legjy,c,jarray,djarray
      dimension ua(0:nmax-1),legjy(0:nmax1-1,0:6,0:3,0:lmaxa1)
      dimension jarray(-lmax1-1:lmax1,0:6),djarray(-lmax1-1:lmax1,0:6)
      external legendrejn

      do n=0,nmax-1
         call legendrejn(ua(n),c,lf(lmaxa),jarray,djarray,lmax1)
         do la=0,lmaxa
            l=max(lf(la),0)
            mult=(-1)^(l+1)
            do a=0,6
               jl(n,a,la)=jarray(l,a)
               yl(n,a,la)=mult*jarray(-l-1,a)
               djl(n,a,la)=djarray(l,a)
               dyl(n,a,la)=mult*djarray(-l-1,a)
            end do
         end do
      end do

      return
      end

@ Figure the workspace for |initjy| (|jarray| and |djarray|).
@<Index spec...@>=
      workalloc(2*(lf(lmaxa)+2+1)*7 * 2)

@ Calculate the potentials and their derivates for the $l$th Legendre
harmonic.

These are some definitions so we don't have to pass a whole slew of arrays
around to subroutines.

@m j(n,a)      legjy(n,a,0)     /* Short-hand expressions for $j$, $y$, their
 derivatives and integrals. */
@m y(n,a)      legjy(n,a,1)
@m dj(n,a)     legjy(n,a,2)
@m dy(n,a)     legjy(n,a,3)
@m ij(n,a)     ijy(n,a,0)
@m iy(n,a)     ijy(n,a,1)

@m psi(n,a)    psid(n,a,0)
@m dpsi(n,a)   psid(n,a,1)


@ Subroutine to calculate the potentials.  This is a straight
implementation of equations (31) and (22).

@<Functions...@>=
      subroutine potential(ua,du,c,fm,chil,nmax,nmax1,legjy,psid,fla,ijy)

      implicit_none
      integer p0,p02,p022,p1,p11
      parameter (p0=0, p02=1, p022=2, p1=3, p11=4)
      integer j0,j1,j2,j02,j11,j22,j022
      parameter (j0=0, j1=1, j2=2, j02=3, j11=4, j22=5, j022=6)

      integer nmax,nmax1,n,a
      real ua,du,c,fm,chil,legjy
      dimension ua(0:nmax-1),fm(0:nmax-1),chil(0:nmax-1),
     $     legjy(0:nmax1-1,0:6,0:3)
      real psid
      dimension psid(0:nmax1-1,0:4,0:1)
      real fla,ijy
      dimension fla(0:nmax-1),ijy(0:nmax1,0:6,0:1)

      @<Define various indefinite integrals@>

      do n=0,nmax-1
         psi(n,p0) = y(n,j0)*ij(n,j0) + iy(n,j0)*j(n,j0)
         psi(n,p02) = y(n,j0)*ij(n,j02) + y(n,j02)*ij(n,j2) +
     $        iy(n,j0)*j(n,j02) + iy(n,j02)*j(n,j2)
         psi(n,p022) = y(n,j0)*ij(n,j022) + y(n,j02)*ij(n,j22) +
     $        y(n,j022)*ij(n,j2) + @|
     $        iy(n,j0)*j(n,j022) + iy(n,j02)*j(n,j22) + iy(n,j022)*j(n,j2)
         psi(n,p1) = y(n,j1)*ij(n,j1) + iy(n,j1)*j(n,j1)
         psi(n,p11) = y(n,j1)*ij(n,j11) + y(n,j11)*ij(n,j1) +
     $        iy(n,j1)*j(n,j11) + iy(n,j11)*j(n,j1)

         dpsi(n,p0) = dy(n,j0)*ij(n,j0) + iy(n,j0)*dj(n,j0)
         dpsi(n,p02) = dy(n,j0)*ij(n,j02) + dy(n,j02)*ij(n,j2) +
     $        iy(n,j0)*dj(n,j02) + iy(n,j02)*dj(n,j2)
         dpsi(n,p022) = dy(n,j0)*ij(n,j022) + dy(n,j02)*ij(n,j22) +
     $        dy(n,j022)*ij(n,j2) + @|
     $        iy(n,j0)*dj(n,j022) + iy(n,j02)*dj(n,j22) + iy(n,j022)*dj(n,j2)
         dpsi(n,p1) = dy(n,j1)*ij(n,j1) + iy(n,j1)*dj(n,j1)
         dpsi(n,p11) = dy(n,j1)*ij(n,j11) + dy(n,j11)*ij(n,j1) +
     $        iy(n,j1)*dj(n,j11) + iy(n,j11)*dj(n,j1)
      end do

      return
      end

@ Workspace for |potential| (|fla| and |ijy|).
@<Index spec...@>=
      workalloc(nmax + (nmax1+1)*7*2)

@ Define various indefinite integrals needed in the definition of the
potentials.  First do the integration to cell edges and then interpolate to
the cell centers.

@<Define various indefinite integrals@>=
      do n=0,nmax-1
         fla(n)=fm(n)*chil(n)*ua(n)^2/sqrt(1+(ua(n)/c)^2)*du/2
      end do

      do a=0,6
         do n=0,nmax-1                /* Define the integrand */
            ij(n+1,a)=fla(n)*j(n,a)
            iy(n,a)=fla(n)*y(n,a)
         end do
         ij(0,a)=0              /* Define the boundary condtions */
         iy(nmax,a)=0
      end do

      do n=1,nmax          /* Integrate to cell edges */
         do a=0,6
            ij(n,a)=ij(n-1,a)+ij(n,a)
         end do
      end do
      do n=nmax-1,0,-1
         do a=0,6
            iy(n,a)=iy(n+1,a)+iy(n,a)
         end do
      end do

      do a=0,6
         do n=0,nmax-1
            ij(n,a)=ij(n,a)+ij(n+1,a)
            iy(n,a)=iy(n,a)+iy(n+1,a)
         end do
      end do

@* Legendre harmonic decomposition.  Define arrays for Legendre harmonics
and decompostion of $\chi$.

@<Glob...@>=
      real chil,legp      /* Legendre decomposition of $\chi$ */
      dimension chil(0:nmax1-1,1:lmaxa1),legp(0:imax1-1,1:lmaxa1)
      external initp

@
@<Index decl...@>=
      integer ptr(chil),ptr(legp)

@
@<Index spec...@>=
      alloc(chil,nmax1*lmaxa1)
      alloc(legp,imax1*lmaxa1)

@ Initialize |legp|.

@<Init...@>=
      call initp(csth0a,imax,imax1,lmaxa,lmaxa1,legp,work)

@ Define |initp|.  Use the recursion formula for Legendre harmonics
$$ (l+1) P_{l+1}(\mu) = (2l+1)\mu P_l(\mu) -  l P_{l-1}(\mu). $$

@<Functions...@>=
      subroutine initp(x,imax,imax1,lmaxa,lmaxa1,legp,leg0)

      implicit_none
      integer imax,imax1,lmaxa,lmaxa1,i,l,la
      real x,legp
      dimension x(0:imax-1),legp(0:imax1-1,1:lmaxa1)
      real leg0
      dimension leg0(0:imax-1)

      if (lmaxa<1) return    /* Nothing to do */

      la=1
      do i=0,imax-1
         leg0(i)=1
         legp(i,la)=x(i)
      end do

      do la=2,lmaxa
         l=lf(la)
         do i=0,imax-1
            leg0(i)=((2*l-3)*x(i)*legp(i,la-1)-(l-2)*leg0(i))/(l-1)
                        /* $P_{l-1}$ */
            legp(i,la)=((2*l-1)*x(i)*leg0(i)-(l-1)*legp(i,la-1))/l
                        /* $P_{l}$ */
         end do
      end do

      return
      end

@ Workspace for |initp| (|leg0|).
@<Index spec...@>=
      workalloc(imax)

@ Decompose $\chi$ into Legendre components.  We have
$$
\chi_l(u)=\frac{2l+1}{2}\int_0^{\pi} \chi(u,\theta_0) P_l(\cos\theta_0) 
\sin\theta_0\, d\theta_0.
$$
Since $\chi$ is odd in $\theta_0$ and it vanishes for the trapped
population this becomes
$$
\chi_l(u)=(2l+1)\int_0^{\theta_{0,\rm max}} \chi(u,\theta_0)  P_l(\cos\theta_0)
\sin\theta_0\, d\theta_0.
$$


@<Functions...@>=
      subroutine decomp(chi,nmax,nmax1,imax,imax1,lmaxa,lmaxa1,
     $    th0a,dth0,legp,chil)

      implicit_none
      integer nmax,nmax1,imax,imax1,lmaxa,lmaxa1,la,l,n,i
      real chi,th0a,dth0,legp,chil,sn
      dimension chi(0:imax1-1,0:nmax-1),th0a(0:imax-1),
     $     legp(0:imax1-1,1:lmaxa1),chil(0:nmax1-1,1:lmaxa1)

      do la=1,lmaxa
         l=lf(la)

         do n=0,nmax-1
            chil(n,la)=0
         end do

         do i=0,imax-1
            sn=sin(th0a(i))*legp(i,la)
            do n=0,nmax-1
               chil(n,la)=chil(n,la)+chi(i,n)*sn
            end do
         end do

         do n=0,nmax-1
            chil(n,la)=(2*l+1)*dth0*chil(n,la)
         end do

      end do

      return
      end

@* Coefficients for isotropic collision operator.  I.e., collisions off a
Maxwellian background.

@<Glob...@>=
      real t,ze,zi              /* Plasma parameters */

      real duug,fu,dtt
      dimension duug(0:nmax),fu(0:nmax-1),dtt(0:nmax-1)

      real fm      /* Maxwellian  */
      dimension fm(0:nmax-1)

      real psid         /* Arrays for potentials and their derivatives */
      dimension psid(0:nmax1-1,0:4,0:1)

      external maxwellian,isotropic

@
@<Index decl...@>=
      integer ptr(duug),ptr(fu),ptr(dtt),ptr(fm),ptr(psid)

@
@<Index spec...@>=
      alloc(duug,nmax+1)
      alloc(fu,nmax)
      alloc(dtt,nmax)
      alloc(fm,nmax)
      alloc(psid,nmax1*5*2)

@ Initialize isotropic terms.

@<Init...@>=
      call maxwellian(ua,du,c,nmax,fm,t)

      do n=0,nmax-1
         chil(n,1)=1
      end do

      call isotropic(ua,du,c,fm,chil(0,1),nmax,nmax1,legjy(0,0,0,0),psid,
     $     ze,zi,duug,fu,dtt,work)

@ This function calculates $D_{uu}$, $D_{\theta\theta}$ and $F_u$ from the
potentials.  This is a transcription of equations (33).

@<Functions...@>=
      subroutine isotropic(ua,du,c,fm,chil,nmax,nmax1,legjy,psid,
     $    ze,zi,duug,fu,dtt,work)

      implicit_none
      integer p0,p02,p022,p1,p11
      parameter (p0=0, p02=1, p022=2, p1=3, p11=4)

      integer nmax,nmax1,n
      real ua,c,psid,duug,fu,dtt,pi,g,u,ze,zi,du,fm,chil,legjy
      dimension ua(0:nmax-1),fm(0:nmax-1),chil(0:nmax-1),
     $     legjy(0:nmax1-1,0:6,0:3),psid(0:nmax1-1,0:4,0:1),
     $     duug(0:nmax),fu(0:nmax-1),dtt(0:nmax-1)
      real work
      dimension work(0:nmax+(nmax1+1)*7*2-1)
      external potential

      call potential(ua,du,c,fm,chil,nmax,nmax1,legjy,psid,work(0),work(nmax))

      pi=4*atan(const(1.0))
      do n=0,nmax-1
         u=ua(n)
         g=sqrt(1+(u/c)^2)
         duug(n)=ze * 4*pi*g/u*(2*g^2*dpsi(n,p02)-u*psi(n,p0)- @|
     $        8*g^2/c^2*dpsi(n,p022)+8*u/c^4*psi(n,p022))
         dtt(n)=ze * 4*pi/(g*u)*(-g^2*dpsi(n,p02)-u/c^2*psi(n,p02)+ @|
     $        4*g^2/c^2*dpsi(n,p022)-4*u/c^4*psi(n,p022))
     $        + zi * g/(2*u)
         fu(n)=ze * 4*pi*g*(-dpsi(n,p1)+2/c^2*dpsi(n,p11))
      end do

 /* Interpolate |duug| to cell edges.  Assume a constant first derivative at
 |umax| and a constant value at $u=0$. */
      duug(nmax)=(3*duug(nmax-1)-duug(nmax-2))/2
      do n=nmax-1,1,-1
         duug(n)=(duug(n)+duug(n-1))/2
      end do
      duug(0)=duug(0)

      return
      end


@ Define a Maxwellian distribution
$$
f_m(u) \propto \exp\biggl(-\frac{mc^2\gamma}{T}\biggr).
$$

@<Functions...@>=
      subroutine maxwellian(ua,du,c,nmax,fm,t)

      implicit_none
      integer nmax,n
      real ua,du,c,fm,t,sum,pi
      dimension ua(0:nmax-1),fm(0:nmax-1)

      sum=0
      do n=0,nmax-1
         fm(n)=exp(-ua(n)^2/(sqrt(1+(ua(n)/c)^2)+1)/t)
         sum=sum+ua(n)^2*fm(n)
      end do
      pi=4*atan(const(1.0))
      sum=4*pi*du*sum
      do n=0,nmax-1
         fm(n)=fm(n)/sum        /* Normalize to unit density */
      end do

      return
      end

@* Reaction terms.  Calculate the term
$$
I(\chi_l(u))=
\frac{C^{e/e}(f_{m}(u),f_{m}(u)\chi_{l}(u)P_l(\cos\theta))}
{f_{m}(u)P_l(\cos\theta)},
$$
using equation (37).  In fact we use a generalization of this equation for
arbitrary $l$.  This term is added to the |reactl|.

@<Functions...@>=
      subroutine reaction(ua,du,c,fm,t,chil,nmax,nmax1,legjy,psid,ze,reactl,l,
     $    work)

      implicit_none
      integer p0,p02,p022,p1,p11
      parameter (p0=0, p02=1, p022=2, p1=3, p11=4)

      integer nmax,nmax1,n,l
      real ua,du,c,legjy,psid,ze,reactl,fm,chil,t,pi,u,g
      dimension ua(0:nmax-1),legjy(0:nmax1-1,0:6,0:3),psid(0:nmax1-1,0:4,0:1),
     $     reactl(0:nmax-1),fm(0:nmax-1),chil(0:nmax-1)
      real work
      dimension work(0:nmax+(nmax1+1)*7*2-1)
      external potential

      call potential(ua,du,c,fm,chil,nmax,nmax1,legjy,psid,work(0),work(nmax))

      pi=4*atan(const(1.0))
      do n=0,nmax-1
         u=ua(n)
         g=sqrt(1+(u/c)^2)
         reactl(n)=reactl(n)+ze*4*pi*(
     $        fm(n)*chil(n)/g- @|
     $        u/t*dpsi(n,p1)-2/(c^2*g)*psi(n,p1)+2*u/(c^2*t)*dpsi(n,p11)+ @|
     $        u/t*dpsi(n,p0)-(u^2/(g*t^2)-1/t)*psi(n,p0)+ @|
     $        (2*g*u/t^2-2*u/(c^2*t))*dpsi(n,p02)-
     $        (l*(l+1)/(g*t^2)-2/(c^2*t))*psi(n,p02)- @|
     $        8*g*u/(c^2*t^2)*dpsi(n,p022)+
     $        4*((l+2)*(l-1)/g+2*g)/(c^2*t^2)*psi(n,p022)
     $        )
      end do

      return
      end

@ Sum up reaction terms
$$ \frac1{\lambda_1}
  \sum_l^{l_{\rm max}}\sum_{l',l''\ \rm odd}^l
  H_{l,l',l''} P_{l',0} I_l(\chi_{l'',0}).
$$
We can rewrite this as
$$ \frac1{\lambda_1}
  \sum_{l'=1}^{l_{\rm max}} P_{l',0}
    \sum_{l=l'}^{l_{\rm max}}
      I_l\biggl( \sum_{l''=1}^l H_{l,l',l''} \chi_{l'',0} \biggr).
$$

@<Functions...@>=
      subroutine reactall(ua,du,c,fm,t,chil,nmax,nmax1,imax,imax1,lam1a,
     $    legjy,psid,legp,harr,ze,lmaxa,lmaxa1,react,reactl,chila,work)

      implicit_none
      integer nmax,nmax1,imax,imax1,lmaxa,lmaxa1
      real ua,du,c,fm,t,chil,lam1a,legjy,psid,legp,harr,ze,react
      dimension ua(0:nmax-1),fm(0:nmax-1),chil(0:nmax1-1,1:lmaxa1),
     $     lam1a(0:imax-1),legjy(0:nmax1-1,0:6,0:3,1:lmaxa1),
     $     psid(0:nmax1-1,0:4,0:1),legp(0:imax1-1,1:lmaxa1),
     $     hv(lmaxa1,lmaxa1,lmaxa1),react(0:imax1-1,0:nmax-1)
      integer la,la1,la2,n,i
      real reactl,chila,work
      dimension reactl(0:nmax-1),chila(0:nmax-1),work(0:nmax+(nmax1+1)*7*2-1)
      external reaction

      do n=0,nmax-1
         do i=0,imax-1
            react(i,n)=0
         end do
      end do

      do la1=1,lmaxa
         @<Compute innermost sums@>
         do n=0,nmax-1
            do i=0,imax-1
               react(i,n)=react(i,n)+reactl(n)*legp(i,la1)
            end do
         end do
      end do

      do n=0,nmax-1
         do i=0,imax-1
            react(i,n)=react(i,n)/lam1a(i)
         end do
      end do

      return
      end

@ Workspace for |reactall| (|reactl| and |chila|) and for |potential|
(|fla| and |ijy|).
@<Index spec...@>=
      workalloc(nmax * 2 + nmax + (nmax1+1)*7*2)

@ Do the innermost two sums.

@<Compute innermost sums@>=
         do n=0,nmax-1
            reactl(n)=0
         end do
         do la=la1,lmaxa
            do n=0,nmax-1
               chila(n)=0
            end do
            do la2=1,la
               do n=0,nmax-1
                  chila(n)=chila(n)+hv(la,la1,la2)*chil(n,la2)
               end do
            end do
            call reaction(ua,du,c,fm,t,chila,nmax,nmax1,
     $           legjy(0,0,0,la),psid,ze,reactl,lf(la),work)
         end do

@* Inhomogeneous magnetic field.

@<Glob...@>=
      real eps,deltab               /* Inverse aspect ratio */
      external deltab           /* Variation of magnetic field */

@ Define magnetic field variation,
$$b=\frac{B(\phi)}{B(0)},\qquad \delta b=b-1,$$
where $\phi$ is the poloidal angle.  We assume that
$$b=\frac{1+\epsilon}{1+\epsilon\cos\phi};\qquad
\delta b=\frac{2\epsilon\sin^2(\phi/2)}{1+\epsilon\cos\phi}.$$

@<Functions...@>=
      function deltab(phi,eps)

      implicit_none
      real deltab,phi,eps

      deltab=2*eps*sin(phi/2)^2/(1+eps*cos(phi))

      return
      end

@ Arrays to hold various averaged quantities.


@<Glob...@>=
      integer kmax
      real lam1,lam2,lam1a,lam2g,harr
      dimension lam1a(0:imax-1),lam2g(0:imax),hv(lmaxa1,lmaxa1,lmaxa1)
      external lam1,lam2,hinit

@
@<Index decl...@>=
      integer ptr(lam1a),ptr(lam2g),ptr(harr)

@
@<Index spec...@>=
      alloc(lam1a,imax)
      alloc(lam2g,imax+1)
      alloc(harr,hind(lmaxa1,lmaxa1,lmaxa1))

@ Initialize these arrays.

@<Init...@>=
      do i=0,imax-1
         lam1a(i)=lam1(eps,th0a(i),kmax)
      end do

      do i=0,imax
         lam2g(i)=lam2(eps,th0g(i),kmax)
      end do

      call hinit(eps,kmax,lmaxa1,harr,work)

@ Define averages over the magnetic field.  We need three:
$$\eqalign{
\lambda_1&=\frac 1{2\pi} \int d\phi \frac{\cos\theta_0}{\cos\theta},\cr
\lambda_2&=\frac 1{2\pi}
\int \frac{d\phi}{b} \frac{\cos\theta}{\cos\theta_0},\cr
H_{l,l',l''}&=\frac 1{2\pi} \int d\phi\,b G_{l,l'} G_{l,l''},\cr
}$$
where
$$\sin^2\theta=b \sin^2\theta_0,$$
and hence
$$\cos\theta=\cos\theta_0\sqrt{1-\delta b\tan^2\theta_0}.$$
The following functions are specialized to the case of passing particles
only.

@ Define $\lambda_1$.

@<Functions...@>=
      function lam1(eps,th0,kmax)

      implicit_none
      integer kmax,k
      real lam1,eps,th0,phi,deltab,deltabv,sum,costh,pi
      external deltab

      pi=4*atan(const(1.0))
      sum=0
      do k=0,kmax-1
         phi=2*pi*(k+const(0.25))/kmax
         deltabv=deltab(phi,eps)
         costh=sqrt(1-min(const(1.0),deltabv*tan(th0)^2))
         sum=sum+1/costh
      end do
      lam1=sum/kmax

      return
      end

@ Define $\lambda_2$.

@<Functions...@>=
      function lam2(eps,th0,kmax)

      implicit_none
      integer kmax,k
      real lam2,eps,th0,phi,deltab,deltabv,sum,costh,pi
      external deltab

      pi=4*atan(const(1.0))
      sum=0
      do k=0,kmax-1
         phi=2*pi*(k+const(0.25))/kmax
         deltabv=deltab(phi,eps)
         costh=sqrt(1-min(const(1.0),deltabv*tan(th0)^2))
         sum=sum+costh/(1+deltabv)
      end do
      lam2=sum/kmax

      return
      end

@ Define $H_{l,l',l''}$ via
$$
  H_{l,l',l''}=\frac{2l+1}{2l''+1}
\lave{ b G_{l,l'} G_{l,l''}}.
$$

|gv| allows us to index into |garr| conveniently.  Enough space is allowed
to accomodate |gv(la,0)|, |gv(la,la+1)|, and |gv(la,la+2)|.  These can be
initialized to $0$ thereby simplifying the recursion for $G$.

$H$ is defined out to |lmaxa1| because we need $H_{1,1,1}$ even when
$l_{\rm max}=0$.

@<Functions...@>=
      subroutine hinit(eps,kmax,lmaxa1,harr,garr)

      implicit_none
      integer kmax,k,lmaxa1,la,la1,la2,l,l1
      real harr,eps,phi,deltab,b,pi,garr
      dimension hv(lmaxa1,lmaxa1,lmaxa1)
      dimension gv(lmaxa1,lmaxa1)
      external deltab

      pi=4*atan(const(1.0))

      @<Set initialial values for |harr| and |garr|@>

      do k=0,kmax-1
         phi=2*pi*(k+const(0.25))/kmax
         b=1+deltab(phi,eps)
         @<Fill |garr|@>
         do la=1,lmaxa1
            do la1=1,la
               do la2=1,la
                  hv(la,la1,la2)=hv(la,la1,la2)+b*gv(la,la1)*gv(la,la2)
               end do
            end do
         end do
      end do

      do la=1,lmaxa1
         do la1=1,la
            do la2=1,la
               hv(la,la1,la2)=(2*lf(la)+1)/areal(2*lf(la2)+1)*
     $              hv(la,la1,la2)/kmax
            end do
         end do
      end do

      return
      end

@ Workspace for |hinit| (|garr|).
@<Index spec...@>=
      workalloc(gind(lmaxa1,lmaxa1))

@ Initialize |harr| and set the edge values of |garr| to zero.

@<Set initialial values for |harr| and |garr|@>=
      do la=1,lmaxa1
         do la1=1,la
            do la2=1,la
               hv(la,la1,la2)=0
            end do
         end do
      end do

      do la=-2,lmaxa1-1
         gv(la,0)=0
      end do
      do la=-1,lmaxa1-1
         gv(la,la+1)=0
         gv(la,la+2)=0
      end do
      gv(1,1)=1

@ Use the recursion relation for $G$ to fill |garr|.

@<Fill |garr|@>=
      do la=2,lmaxa1
         l=lf(la)
         do la1=1,la
            l1=lf(la1)
            gv(la,la1)=(2*(2*l-3)*(l^2-3*l+1)*gv(la-1,la1) - @|
     $           (l-3)*(l-2)*(2*l-1)*gv(la-2,la1)) / (l*(l-1)*(2*l-5)) + @|
     $           (2*l-3)*(2*l-1)*b*
     $           ( (l1+1)*(l1+2) * gv(la-1,la1+1)/ ((2*l1+3)*(2*l1+5)) - @|
     $           2*(l1^2+l1-1) * gv(la-1,la1)/((2*l1-1)*(2*l1+3)) + @|
     $           l1*(l1-1) * gv(la-1,la1-1)/((2*l1-3)*(2*l1-1)) )/(l*(l-1))
         end do
      end do

@* Residue calculation.  The array |resid| is already initialized with the
reaction term.  Two methods for the calculation are given here: the
straightforward coding of the residue term, which is slow but easy to
check; and a faster method using the arrays |uinv| and |th0inv| which are
used in the tridiagonal elimination.

@m RESIDUE 0 /* Should we calculate the residue in the straightforward way? */

@<Functions...@>=
@#if RESIDUE
      subroutine residue(du,ua,ug,nmax,dth0,th0a,th0g,csth0a,imax,imax1,
     $    chi,duug,fu,dtt,c,lam1a,lam2g,rho,resid)

      implicit_none
      integer nmax,imax,imax1,n,i
      real du,ua,ug,dth0,th0a,th0g,csth0a,chi,duug,fu,dtt,c,
     $     lam1a,lam2g,rho,resid
      dimension ua(0:nmax-1),ug(0:nmax),th0a(0:imax-1),th0g(0:imax),
     $     csth0a(0:imax-1),
     $     chi(0:imax1-1,0:nmax-1),duug(0:nmax),fu(0:nmax-1),dtt(0:nmax-1),
     $     lam1a(0:imax-1),lam2g(0:imax),resid(0:imax1-1,0:nmax-1)

      @<Energy diffusion term@>

      @<Friction term@>

      @<Pitch-angle scattering term@>

      @<Electric field term@>

 /* Normalize to pitch-angle diffusion coefficient. */
      do n=0,nmax-1
         do i=0,imax-1
            resid(i,n)=resid(i,n)*ua(n)^2/dtt(n)
         end do
      end do

      return
      end
@#endif

@ The faster way.

@<Functions...@>=
@#if !RESIDUE
@#if VAX
      options /check=nobounds
@#endif

      subroutine residue(ua,nmax,nmax1,csth0a,imax,imax1,chi,dtt,c,lam1a,
     $    uinv,th0inv,resid)

      implicit_none
      integer nmax,nmax1,imax,imax1,n,i
      real ua,csth0a,chi,dtt,c,lam1a,resid,uinv,th0inv
      dimension ua(0:nmax-1),chi(0:imax1-1,0:nmax-1),csth0a(0:imax-1),
     $     dtt(0:nmax-1),lam1a(0:imax-1),resid(0:imax1-1,0:nmax-1),
     $     uinv(0:nmax1-1,-1:1),th0inv(0:imax1-1,-1:1)

      @<Electric field term@>

 /* Normalize to pitch-angle diffusion coefficient and add in
    differential terms. */
      do n=0,nmax-1
         do i=0,imax-1
            resid(i,n)=resid(i,n)*ua(n)^2/dtt(n)+ @|
     $           chi(i,n-1)*uinv(n,-1)+chi(i,n)*uinv(n,0)+
     $           chi(i,n+1)*uinv(n,1)+ @|
     $           chi(i-1,n)*th0inv(i,-1)+chi(i,n)*th0inv(i,0)+
     $           chi(i+1,n)*th0inv(i,1)
         end do
      end do

      return
      end
@#endif

@ Calculate contributions to the residue from the energy diffusion term
$$\frac{1}{u^2}\frac{\d}{\d u} u^2 D_{uu} \frac{\d \chi}{\d u}.$$
Boundary at $u=0$ takes care of itself.  Boundary at $u=u_{\rm max}$
handled by specifying
|chi(i,nmax)=chi(i,nmax-1)+rho*(chi(i,nmax-1)-chi(i,nmax-2))|.  For
$\rho=1$, we have $\chi''(u_{\rm max})=0$, and for $\rho=0$, we have
$\chi'(u_{\rm max})=0$.

@<Energy diffusion term@>=
      do i=0,imax-1
         n=0
         resid(i,n)=resid(i,n) + @|
     $        (ug(n+1)^2*duug(n+1)*(chi(i,n+1)-chi(i,n))/du)/(ua(n)^2*du)
         do n=1,nmax-2
            resid(i,n)=resid(i,n) + @|
     $           (ug(n+1)^2*duug(n+1)*(chi(i,n+1)-chi(i,n))/du - @|
     $           ug(n)^2*duug(n)*(chi(i,n)-chi(i,n-1))/du)/(ua(n)^2*du)
         end do
         n=nmax-1
         resid(i,n)=resid(i,n) + @|
     $        (ug(n+1)^2*duug(n+1)*rho*(chi(i,n)-chi(i,n-1))/du - @|
     $        ug(n)^2*duug(n)*(chi(i,n)-chi(i,n-1))/du)/(ua(n)^2*du)
      end do

@ Calculate contributions to the residue from the friction term
$$F_{u} \frac{\d \chi}{\d u}.$$
Boundary at $u=0$ handled by assuming $\chi$ is odd in $u$, i.e.,
|chi(i,-1)=-chi(i,0)|.  Boundary at $u=u_{\rm max}$ handled in the same way
as in the calculation of the energy diffusion term.

@<Friction term@>=
      do i=0,imax-1
         n=0
         resid(i,n)=resid(i,n)+fu(n)*(chi(i,n+1)+chi(i,n))/(2*du)
         do n=1,nmax-2
            resid(i,n)=resid(i,n)+fu(n)*(chi(i,n+1)-chi(i,n-1))/(2*du)
         end do
         n=nmax-1
         resid(i,n)=resid(i,n)+fu(n)*(1+rho)*(chi(i,n)-chi(i,n-1))/(2*du)
      end do

@ Calculate contributions to the residue from the pitch-angle scattering
term
$$\frac{D_{\theta\theta}}{u^2}
\frac{1}{\lambda_1\sin\theta_0}\frac{\d}{\d\theta_0} \sin\theta_0 \lambda_2
\frac{\d}{\d\theta_0}.$$
Boundary at $\theta_0=0$ takes care of itself.  At $\theta_{0,\rm max}
=\theta_t$, assume that $\chi=0$.

@<Pitch-angle scattering term@>=
      do n=0,nmax-1
         i=0
         resid(i,n)=resid(i,n)+dtt(n)/ua(n)^2 * @|
     $        (sin(th0g(i+1))*lam2g(i+1)*(chi(i+1,n)-chi(i,n))/dth0)/
     $        (lam1a(i)*sin(th0a(i))*dth0)
         do i=1,imax-2
            resid(i,n)=resid(i,n)+dtt(n)/ua(n)^2 * @|
     $           (sin(th0g(i+1))*lam2g(i+1)*(chi(i+1,n)-chi(i,n))/dth0 - @|
     $           sin(th0g(i))*lam2g(i)*(chi(i,n)-chi(i-1,n))/dth0)/
     $           (lam1a(i)*sin(th0a(i))*dth0)
         end do
         i=imax-1
         resid(i,n)=resid(i,n)+dtt(n)/ua(n)^2 * @|
     $        (sin(th0g(i+1))*lam2g(i+1)*(-chi(i,n))/(dth0/2) - @|
     $        sin(th0g(i))*lam2g(i)*(chi(i,n)-chi(i-1,n))/dth0)/
     $        (lam1a(i)*sin(th0a(i))*dth0)
      end do

@ Finally add in the electric field term:
$$\frac{v \cos\theta_0}{\lambda_1}. $$

@<Electric field term@>=
      do n=0,nmax-1
         do i=0,imax-1
            resid(i,n)=resid(i,n)+
     $           csth0a(i)/lam1a(i)*ua(n)/sqrt(1+(ua(n)/c)^2)
         end do
      end do

@ Calculate norm of residue.

@<Functions...@>=
      subroutine error(du,ua,umax,nmax,dth0,th0a,th0max,imax,imax1,resid,chi,
     $    abserr,relerr,aerr,rerr)

      implicit_none
      integer nmax,imax,imax1,i,n
      real du,ua,umax,dth0,th0a,th0max,resid,chi,abserr,relerr
      dimension ua(0:nmax-1),th0a(0:imax-1),resid(0:imax1-1,0:nmax-1),
     $     chi(0:imax1-1,0:nmax-1)
      real aerr,rerr
      dimension aerr(0:imax-1),rerr(0:imax-1)

      do i=0,imax-1
         aerr(i)=0
         rerr(i)=0
      end do

      do n=0,nmax-1
         do i=0,imax-1
            aerr(i)=aerr(i)+resid(i,n)^2*ua(n)^2
            rerr(i)=rerr(i)+(resid(i,n)/max(epsilon,abs(chi(i,n))))^2*ua(n)^2
         end do
      end do

      abserr=0
      relerr=0
      do i=0,imax-1
         abserr=abserr+sin(th0a(i))*aerr(i)
         relerr=relerr+sin(th0a(i))*rerr(i)
      end do

 /* Normalize to total volume of velocity space */
      abserr=sqrt(abserr*du*dth0/(umax^3/3*(1-cos(th0max))))
      relerr=sqrt(relerr*du*dth0/(umax^3/3*(1-cos(th0max))))

      return
      end

@ Workspace for |error| (|aerr| and |rerr|).
@<Index spec...@>=
      workalloc(2*imax)

@* Matrices for inversion.

@<Glob...@>=
      real uinv,th0inv,rho
      dimension uinv(0:nmax1-1,-1:1),th0inv(0:imax1-1,-1:1)
      external matrixinit

@
@<Index decl...@>=
      integer ptr(uinv),ptr(th0inv)

@
@<Index spec...@>=
      alloc(uinv,nmax1*3)
      alloc(th0inv,imax1*3)

@ Initialize matrices.

@<Init...@>=
      call matrixinit(du,ua,ug,nmax,nmax1,dth0,th0a,th0g,imax,imax1,
     $    duug,fu,dtt,lam1a,lam2g,uinv,th0inv,rho)

@ Calculate inversion matrices.  These have to be compatible with
subroutine |residue|.

@<Functions...@>=
      subroutine matrixinit(du,ua,ug,nmax,nmax1,dth0,th0a,th0g,imax,imax1,
     $    duug,fu,dtt,lam1a,lam2g,uinv,th0inv,rho)

      implicit_none
      integer nmax,nmax1,imax,imax1,n,i
      real du,ua,ug,dth0,th0a,th0g,duug,fu,dtt,lam1a,lam2g,
     $     uinv,th0inv,rho
      dimension ua(0:nmax-1),ug(0:nmax),th0a(0:imax-1),th0g(0:imax),
     $     duug(0:nmax),fu(0:nmax-1),dtt(0:nmax-1),
     $     lam1a(0:imax-1),lam2g(0:imax),
     $     uinv(0:nmax1-1,-1:1),th0inv(0:imax1-1,-1:1)

      @<Matrices for $u$ inversion@>

      @<Matrices for $\theta_0$ inversion@>

      return
      end

@ Define matrices for inversion of diffusion and friction operator.  These
are one-dimensional because the operator is independent of $\theta_0$.

@<Matrices for $u$ inversion@>=
      do n=0,nmax-1
         uinv(n,-1)=ug(n)^2*duug(n)/(ua(n)^2*du^2)
         uinv(n,1)=ug(n+1)^2*duug(n+1)/(ua(n)^2*du^2)
         uinv(n,0)=-uinv(n,-1)-uinv(n,1)
         uinv(n,-1)=uinv(n,-1)-fu(n)/(2*du)
         uinv(n,1)=uinv(n,1)+fu(n)/(2*du)
      end do
      n=0                       /* Odd at $u=0$ */
      uinv(n,0)=uinv(n,0)-uinv(n,-1)
      uinv(n,-1)=0
      n=nmax-1                  /* Apply boundary condition at $u_{\rm max}$ */
      uinv(n,-1)=uinv(n,-1)-rho*uinv(n,1)
      uinv(n,0)=uinv(n,0)+(1+rho)*uinv(n,1)
      uinv(n,1)=0

      do n=0,nmax-1
         uinv(n,-1)=uinv(n,-1)*ua(n)^2/dtt(n)
         uinv(n,0)=uinv(n,0)*ua(n)^2/dtt(n)
         uinv(n,1)=uinv(n,1)*ua(n)^2/dtt(n)
      end do

@ Define matrices for inversion of pitch-angle scattering operator.  These
are one-dimensional because, after divided out $D_{\theta\theta,0}/u^2$,
the operator is independent of $u$.

@<Matrices for $\theta_0$ inversion@>=
      do i=0,imax-1
         th0inv(i,-1)=sin(th0g(i))*lam2g(i)/(lam1a(i)*sin(th0a(i))*dth0^2)
         th0inv(i,1)=sin(th0g(i+1))*lam2g(i+1)/(lam1a(i)*sin(th0a(i))*dth0^2)
         th0inv(i,0)=-th0inv(i,-1)-th0inv(i,1)
      end do
      i=0                       /* Nothing special needs to be done */
      i=imax-1                  /* Zero on trapped passing boundary */
      th0inv(i,0)=th0inv(i,0)-th0inv(i,1)
      th0inv(i,1)=0

@ Solve the tridiagonal equations.

This subroutine solves the tridiagonal equation
$$
x_i-\alpha\Delta t(a_i x_{i-1}+b_i x_i+c_i x_{i+1})=z_i,
$$
for $0\le i< n$, with $a_0=c_{n-1}=0$.
Vectorized to handle $m$ parallel systems.  |ms| is the skipping distance
in the vectorizing direction.  |ns| is the skipping distance in the
solution direction between $x_i$ and $x_{i+1}$.  In this version, |a|, |b|, and
|c| are the same for each system.

Avoid using any working storage by clobbering |z| and using |x| as a
temporary.

@<Functions...@>=
@#if VAX
      options /check=nobounds
@#endif

      subroutine solve(x,ns,n,ms,m,a,b,c,z,dt,alpha)

      implicit_none
      integer n,ns,m,ms,i,k
      real x,a,b,c,z,dt,alpha,dt2,den
      dimension x(0:ms-1,0:*),z(0:ms-1,0:*),a(0:n-1),b(0:n-1),c(0:n-1)

      dt2=-alpha*dt

      i=0
      do k=0,m-1
         den=1+dt2*b(i)
         x(ns*i,k)=z(ns*i,k)/den
         z(ns*i,k)=-dt2*c(i)/den
      end do

      do i=1,n-2
         do k=0,m-1
            den=1+dt2*(b(i)+a(i)*z(ns*(i-1),k))
            x(ns*i,k)=(z(ns*i,k)-dt2*a(i)*x(ns*(i-1),k))/den
            z(ns*i,k)=-dt2*c(i)/den
         end do
      end do

      i=n-1
      do k=0,m-1
         den=1+dt2*(b(i)+a(i)*z(ns*(i-1),k))
         x(ns*i,k)=(z(ns*i,k)-dt2*a(i)*x(ns*(i-1),k))/den
         z(ns*i,k)=0
      end do


 /* The equations have now been converted to the form $x_i=r_i x_{i+1}+s_i$
 (with $r=y$, $s=x$). */

      do i=n-2,0,-1
         do  k=0,m-1
            x(ns*i,k)=z(ns*i,k)*x(ns*(i+1),k)+x(ns*i,k)
         end do
      end do

      return
      end

@* Time advancement.

@<Glob...@>=
      integer tmax,time,tskip
      real resid,abserr,relerr,alpha,dt,resida
      dimension resid(0:imax1-1,0:nmax-1),resida(0:imax1-1,0:nmax-1)
      external decomp,reactall,residue,error,solve

@
@<Index decl...@>=
      integer ptr(resid),ptr(resida)

@
@<Index spec...@>=
      alloc(resid,imax1*nmax)
      alloc(resida,imax1*nmax)

@ Calculate reaction term.

@<Make a step@>=
      call decomp(chi,nmax,nmax1,imax,imax1,lmaxa,lmaxa1,th0a,dth0,legp,chil)
      call reactall(ua,du,c,fm,t,chil,nmax,nmax1,imax,imax1,lam1a,
     $     legjy(0,0,0,1),psid,legp,harr,ze,lmaxa,lmaxa1,resid,
     $     work(0),work(nmax),work(2*nmax))

@ Calculate residue and its norm.

@<Make a step@>=
@#if RESIDUE
      call residue(du,ua,ug,nmax,dth0,th0a,th0g,csth0a,imax,imax1,chi,
     $    duug,fu,dtt,c,lam1a,lam2g,rho,resid)
@#else
      call residue(ua,nmax,nmax1,csth0a,imax,imax1,chi,dtt,c,lam1a,
     $    uinv,th0inv,resid)
@#endif
      if (mod(time,tskip)==0) then
         call error(du,ua,umax,nmax,dth0,th0a,th0max,imax,imax1,resid,chi,
     $        abserr,relerr,work(0),work(imax))
      end if

@ Advance $\chi$.

@<Make a step@>=
      call solve(resida,imax1,nmax,1,imax,uinv(0,-1),uinv(0,0),uinv(0,1),
     $    resid,dt,alpha)
      call solve(resid,1,imax,imax1,nmax,th0inv(0,-1),th0inv(0,0),th0inv(0,1),
     $    resida,dt,alpha)

      do n=0,nmax-1
         do i=0,imax-1
            chi(i,n)=chi(i,n)+dt*resid(i,n)
         end do
      end do

@* Current and conductivity.

@<Glob...@>=
      real curv,current,cond,e
      external current

@ Compute current and conductivity.  The conductivity is define as
$J_{0\parallel}/E_{0\parallel}$ where
$$E_{0\parallel} = B_0\frac{\int E_{\parallel} dl}{\int B\,dl}.$$
Substituting $E_\parallel = (B_\zeta/B) \Phi/(2\pi R)$ and $\Phi=TQ/L$ gives
$E_\parallel = T/\int b \,dl/L$.

@<Make a step@>=
      if (mod(time,tskip)==0) then
         curv=current(du,ua,c,nmax,dth0,th0a,imax,imax1,fm,chi,work)
         e=t/hv(1,1,1)
         cond=curv*zi/(t*sqrt(t)*e)
         write(outunit,*) time,cond,abserr,relerr
      end if

@ Function to calculate current.

@<Functions...@>=
      function current(du,ua,c,nmax,dth0,th0a,imax,imax1,fm,chi,cur)

      implicit_none
      integer nmax,imax,imax1,n,i
      real current,du,ua,c,dth0,th0a,fm,chi
      dimension ua(0:nmax-1),th0a(0:imax-1),
     $     fm(0:nmax-1),chi(0:imax1-1,0:nmax-1)
      real cur,pi
      dimension cur(0:imax-1)

      pi=4*atan(const(1.0))
      do i=0,imax-1
         cur(i)=0
      end do

      do n=0,nmax-1
         do i=0,imax-1
            cur(i)=cur(i)+chi(i,n)*fm(n)*ua(n)**3/sqrt(1+(ua(n)/c)^2)
         end do
      end do

      current=0
      do i=0,imax-1
         current=current+sin(th0a(i))*cos(th0a(i))*cur(i)
      end do

      current=4*pi*du*dth0*current

      return
      end

@ Workspace for |current| (|cur|).
@<Index spec...@>=
      workalloc(imax)

@* Write out results to disk.

@<Write out results@>=
      write(diskunit) nmax,imax,kmax,lmax
      write(diskunit) single(t),single(c2),single(ze),single(zi),single(eps)
      write(diskunit) single(umax),single(th0max),single(du),single(dth0),
     $     single(rho)
      write(diskunit) single(dt),tmax,single(alpha)
      write(diskunit) single(cond),single(abserr),single(relerr)
      write(diskunit) (single(ua(n)),n=0,nmax-1)
      write(diskunit) (single(ug(n)),n=0,nmax)
      write(diskunit) (single(th0a(n)),n=0,imax-1)
      write(diskunit) (single(th0g(n)),n=0,imax)
      write(diskunit) (single(fm(n)),n=0,nmax-1)
      write(diskunit) ((single(chi(i,n)),i=0,imax-1),n=0,nmax-1)

@* Legendre functions.  It is convenient to express the potentials in terms
of a set of functions defined by
$$
\eqalignno{
j\subb l1a(z) &=\sqrt{\frac{\pi}{2z}}P^{-l-1/2}_{a-1/2}(\sqrt{1+z^2}),\cr
y\subb l1a(z) &=(-1)^{l+1}\sqrt{\frac{\pi}{2z}}P^{l+1/2}_{a-1/2}(\sqrt{1+z^2}),
\cr}
$$
where $P^\mu_\nu$ is the associated Legendre function of the first kind.
Multiple index functions are defined recursively by
$$
j\subb l{k+2}{\s aa'} =
\cases{
  \displaystyle\frac{j\subb l{k+1}{\s a}-j\subb l{k+1}{\s a'}}
    {a^2-a'^2},&for $a\ne a'$,\cr
  \displaystyle\pd{j\subb l{k+1}{\s a}}{(a^2)},
    &for $a=a'$,\cr
}$$
with $y\subb lk\s$ defined in a similar fashion.
The two families of functions are related by
$$y\subb lk{\s}=(-1)^{l+1}j\subb {-l-1}k{\s}.$$

In fact we shall work with two derived functions, scaled Legendre functions
$\ja{l,a}(z)$ defined by
$$j\subb l1a(z) = \frac{z^l}{(2l+1)!!} \ja{l,a}(z),$$
and normalized Legendre functions defined by
$$ \jn lk{\s}(u) = c^{2k+l-2} j\subb lk{\s}(u/c). $$
For $k=1$, we have
$$\jn\subb l1a(z) = \frac{u^l}{(2l+1)!!} \ja{l,a}(u/c).$$

The methods used for computing these functions are given in the appendix.

@ Scaled Legendre functions.  This subroutine returns the scaled Legendre
functions $\hbox{|jarray(l,a)|}=\ja{l,a}(z)$ for $-|lmax|-1 \le l \le
|lmax|$ and $0\le a \le |amax|$.

@<Functions...@>=
      subroutine legendrejs(z,lmax,amax,jarray,lmax1)

      implicit_none
      integer l,a,lmax,amax,lmax1,lstart
      real z,jarray,z2,g,ja0,ja1,ja2,ratio,sigma
      dimension jarray(-lmax1-1:lmax1,0:amax)

      if (amax<0 || lmax<0) return /* Nothing to do! */
      if (lmax>lmax1) return /* Error: not enough room! */

      z2=z^2
      g=sqrt(1+z2)

      @< Do the cases $-a\le l<a$ @>

      @< Do the cases $l<-a$ @>

      @< Do the cases $l\ge0$, $a=0$ @>

      @< Do the other cases, i.e., $a>0$, $l\ge a$ @>

      return
      end

@ Calculate the scaled Legendre functions for the easy case, namely $-a\le
l<a$.

@< Do the cases $-a\le l<a$ @>=

      do a=1,amax               /* $a=0$ is excluded */
         ja2=1                  /* $\ja{a-1,a}$ */
         ja1=g                  /* $\ja{a-2,a}$ */
         do l=a-3,lmax-1,-1 @/
 /* Get things started by iterating down to $|lmax|-1$ if need be.
 This is from equation (A29) */
            ja0=g*ja1-z2*(l-a+2)*(l+a+2)*ja2/((2*l+3)*(2*l+5))
            ja2=ja1             /* $|ja2| = \ja{l+1,a}$ */
            ja1=ja0             /* $|ja1| = \ja{l,a}$ */
         end do
         l=min(lmax,a-1)
         jarray(l,a)=ja2
         l=l-1
         jarray(l,a)=ja1
         do l=min(lmax,a-1)-2,max(-a,-lmax-1),-1 @/
 /* Main iteration loop---again from equation (A29) */
            jarray(l,a)=g*jarray(l+1,a)-
     $           z2*(l-a+2)*(l+a+2)*jarray(l+2,a)/((2*l+3)*(2*l+5))
         end do
      end do

@ Calculate the scaled Legendre functions for the other easy case, namely
$l<-a$.

@< Do the cases $l<-a$ @>=

      do a=0,min(amax,lmax)
         l=-a-1
         jarray(l,a)=1
         l=l-1
         if (l>=-lmax-1) then
            jarray(l,a)=g
         end if
         do l=-a-3,-lmax-1,-1 @/
 /* Main iteration loop---again from equation (A29) */
            jarray(l,a)=g*jarray(l+1,a)-
     $           z2*(l-a+2)*(l+a+2)*jarray(l+2,a)/((2*l+3)*(2*l+5))
         end do
      end do

@ Calculate $\ja{l\ge0,0}$.  Depending on the magnitude of the argument we
either use backwards recursion (for $z \le |lmax|+1$) or forwards
recursion using $\ja{-1,0}=1$ and $\ja{0,0} = (\sinh^{-1} z)/z$.

@< Do the cases $l\ge0$, $a=0$ @>=
      a = 0
      if (abs(z) <= lmax+1) then @/
 /* Estimate where to start the recursion.  The 10 is a extra guard. */
         lstart = lmax + log(epsilon)/(2*log((epsilon+abs(z))/(g+1))) + 10 @/
 /* Initial guess for ratio $\ja{L'+1,a}/\ja{L',a}$. */
         ratio = 2/(g+1)
         do l=lstart,lmax-1,-1
            ratio=1/(g - z2*(l-a+2)*(l+a+2)*ratio/((2*l+3)*(2*l+5))) @/
 /* At this point, |ratio| is approximation for $\ja{l+1,a}/\ja{l,a}$. */
         end do
         jarray(lmax,a)=ratio
         ja0=1
         do l=lmax-2,-1,-1
            jarray(l+1,a)=ja0
            ja0=g*jarray(l+1,a)-
     $           z2*(l-a+2)*(l+a+2)*jarray(l+2,a)/((2*l+3)*(2*l+5))
         end do
         ja0=1/ja0
         do l=0,lmax
            jarray(l,a)=ja0*jarray(l,a)
         end do
      else
         /* Use forward recurrence for $z>|lmax|+1$. */
         sigma=log(g+abs(z))        /* $\sinh^{-1}(\abs z)$ */
         l=0
         jarray(l,a)=sigma/abs(z)   /* Starting value.  We assume
             that $\ja{-1,0}$ has already been set to $1$. */
         do l=1,lmax @/
 /* Equation (A29) */
            jarray(l,a)=(g*jarray(l-1,a)-jarray(l-2,a))*
     $           (2*l-1)*(2*l+1)/(z2*(l-a)*(l+a))
         end do
      end if

@ Do the other cases, i.e., $a>0$, $l\ge a$.

@< Do the other cases, i.e., $a>0$, $l\ge a$ @>=
      do a=1,min(lmax,amax)
         l=lmax @/
 /* Initialize backwards recursion using equation (A30). */
         jarray(l,a)=((2*l+1)*jarray(l-1,a-1)-
     $    (l+1-a)*g*jarray(l,a-1))/(l+a)
         l=lmax-1
         if (l==0) then         /* Implies that $a=1$.  Use $\ja{-1,0} = 1$ */
            ja0=((2*l+1)-(l+1-a)*g*jarray(l,a-1))/(l+a)
         else
            ja0=((2*l+1)*jarray(l-1,a-1)-
     $           (l+1-a)*g*jarray(l,a-1))/(l+a)
         end if
         do l=lmax-2,a-1,-1
            jarray(l+1,a)=ja0
            ja0=g*jarray(l+1,a)-
     $           z2*(l-a+2)*(l+a+2)*jarray(l+2,a)/((2*l+3)*(2*l+5))
         end do
         ja0=1/ja0
         do l=a,lmax
            jarray(l,a)=ja0*jarray(l,a)
         end do
      end do

@ Legendre functions of the first kind.  This returns $
|jarray(l,a)|=j\subb l1a(z)$ for $0\le l \le |lmax|$ and $0\le a \le
|amax|$.  In fact, we don't need this function; so it's commented out.

@<Functions...@>=
@#if 0
      subroutine legendrej(z,lmax,amax,jarray,lmax1)

      implicit_none
      integer l,a,lmax,amax,lmax1
      real z,jarray,scale
      dimension jarray(-lmax1-1:lmax1,0:amax)
      external legendrejs

      call legendrejs(z,lmax,amax,jarray,lmax1)

      scale=1
      do l=1,lmax
         scale=scale*z/(2*l+1)
         do a=0,amax
            jarray(l,a)=scale*jarray(l,a)
         end do
      end do

      scale=1
      do l=-1,-lmax-1,-1
         scale=scale/z*(2*l+3)
         do a=0,amax
            jarray(l,a)=scale*jarray(l,a)
         end do
      end do

      return
      end
@#endif

@ Calculate normalized Legendre functions needed in computation of
collision operator.  The single-index functions
$\jn l1a(u)$ are given in terms of $\ja{l,a}(u/c)$.
The multiple-index function are given in terms of the single-index
functions by by equations (A26):
$$\eqalign{
\jn l2{02}(u)&=\frac u2 \jn{l+1}11(u),\cr
\jn l2{11}(u)&=\frac u2 \jn{l+1}10(u),\cr
\jn l2{22}(u)&=\frac u2 \jn{l+1}11(u)+\frac{u^2}{2c^2} \jn{l+2}10(u),\cr
\jn l3{022}(u)&=\frac {u^2}8 \jn{l+2}10(u).\cr
}$$
The derivatives are given by equation (A21):
$$
\frac d{du}\jn lk*(u)=
  \frac1{\gamma}\jn{l-1}k*(u)-\frac{l+1}{u}\jn lk*(u).
$$
On output, this returns
$$\eqalign{
\hbox{|jarray(l,0)|} &= \jn l10(u), \cr
\hbox{|jarray(l,1)|} &= \jn l11(u), \cr
\hbox{|jarray(l,2)|} &= \jn l12(u), \cr
\hbox{|jarray(l,3)|} &= \jn l2{02}(u), \cr
\hbox{|jarray(l,4)|} &= \jn l2{11}(u), \cr
\hbox{|jarray(l,5)|} &= \jn l2{22}(u), \cr
\hbox{|jarray(l,6)|} &= \jn l3{022}(u), \cr
\hbox{|djarray(l,a)|} &=\frac d{du}\hbox{|jarray(l,a)|}, \cr
}$$
for $-|lmax|-1 \le l \le |lmax|$ and $0\le a \le |amax|$.  The dimensioning
of the arrays is controlled by |lmax1| which should be at least $|lmax|+2$,
in order to provide additional work space.

@ Here is |legendrejn|.

@<Functions...@>=
      subroutine legendrejn(u,c,lmax,jarray,djarray,lmax1)

      implicit_none
      integer j0,j1,j2,j02,j11,j22,j022
      parameter (j0=0, j1=1, j2=2, j02=3, j11=4, j22=5, j022=6)

      integer l,a,lmax,amax,lmax1
      real u,c,jarray,djarray,scale,g
      dimension jarray(-lmax1-1:lmax1,0:6),djarray(-lmax1-1:lmax1,0:6)
      external legendrejs

      amax=2
      call legendrejs(u/c,lmax+2,amax,jarray,lmax1)

      scale=1
      do l=1,lmax+2
         scale=scale*u/(2*l+1)
         do a=0,amax
            jarray(l,a)=scale*jarray(l,a)
         end do
      end do

      scale=1
      do l=-1,-lmax-2,-1
         scale=scale/u*(2*l+3)
         do a=0,amax
            jarray(l,a)=scale*jarray(l,a)
         end do
      end do

      do l=-lmax-2,lmax         /* Calculate multiple index solutions */
         jarray(l,j02) = u/2*jarray(l+1,1)
         jarray(l,j11) = u/2*jarray(l+1,0)
         jarray(l,j22) = u/2*jarray(l+1,1) + (u/c)^2/2*jarray(l+2,0)
         jarray(l,j022) = u^2/8*jarray(l+2,0)
      end do

      g = sqrt(1 + (u/c)**2)         /* Calculate derivatives */
      do a=0,6
         do l=-lmax-1,lmax
            djarray(l,a) = jarray(l-1,a)/g - (l+1)/u*jarray(l,a)
         end do
      end do

      return
      end

@* Index.
