c ======================================================================
c     usage: partdemo1 <np> [<a>]
c
c     Demonstrates use of 'gft_out_part_xyzf' to output particle data 
c     with a number of scalar functions (two in this case) defined 
c     on the particles. 
c
c     Specifically, this program generates <np> particles randomly 
c     distributed on the surface of a torus, where <a> defines the 
c     "aspect ratio" of the torus, so that as a -> 0, the torus 
c     becomes a ring with unit radius, centred on the origin, and 
c     lying in the z=0 plane.
c
c     The scalar functions defined on the particles are 
c 
c     1) min(x_j,y_j,z_j)
c     2) max(x_j,y_j,z_j)
c
c     where (x_j,y_j,z_j) are the coordinates of the j-th point.
c
c     In addition to 'gft_out_part_xyzf', defined in the 'bbhutil'
c     library, this program references the external routines 'i4arg',
c     'r8arg', 'srand' and 'drand48', which are defined in the 'p410f' 
c     library.
c ======================================================================
      program      partdemo1

      implicit     none

      character*9  cdnm
      parameter  ( cdnm = 'partdemo1' )

      real*8       r8arg
      integer      iargc,  i4arg

      integer      np
      real*8       a  

      integer      npmax
      parameter  ( npmax = 10 000 000 )

      real*8       r(npmax,6)
      real*8       t

      if( iargc() .lt. 1  .or.  iargc() .gt. 2  ) go to 900

      np = i4arg(1,-1)
      if( np .eq. -1 ) go to 900
      if( np .lt. 1  .or.  np .gt. npmax ) then
         write(0,*) 'cdnm: np=', np
         write(0,*) 'cdnm: Exceeds range 1 to ', npmax
         go to 900
      end if
      a = r8arg(2,0.1d0)

      write(0,*) cdnm//': a=', a
      call genpart(r,np,a)
      t = 0.0d0 

      call gft_out_part_xyzf('partdemo1',t,r,np,2)

      stop

 900  continue
         write(0,*) 'usage: '//cdnm//' <np> [<a>]'
      stop
      
      end

c ======================================================================
c     Generates particle coordinates and particle functions, and returns 
c     them in the 2-d array, 'r'.
c ======================================================================
      subroutine genpart(r,np,a)
         implicit     none

         real*8       drand48

         integer      np
         real*8       r(np,5)
         real*8       a

         integer      j
         real*8       pi, th1, th2, xj, yj, zj, nrm

         call srand(np)

         pi = 4.0d0 * atan(1.0d0)
         do j = 1 , np 
            th1 = pi * (2.0d0 * drand48() - 1.0d0)
            th2 = pi * (2.0d0 * drand48() - 1.0d0)
            xj = cos(th1) * (1.0d0 + a * cos(th2))
            yj = sin(th1) * (1.0d0 + a * cos(th2))
            zj = a * sin(th2)
            r(j,1) = xj
            r(j,2) = yj
            r(j,3) = zj
            r(j,4) = max(xj,yj,zj)
            r(j,5) = min(xj,yj,zj)
         end do

         return

      end