c---------------------------------------------------------------------
c     Brute force evolution of charges on sphere.
c
c     Requires input file
c
c        chgsph.control
c---------------------------------------------------------------------
c
      real*8         dsqrt, dvvdot, drand48
      integer        mymod
c
      integer        indlnb,          iargc
      integer        narg,            iarg,        larg
      character*80   argi,            pgname
c
      integer        mxspos,          mxch 
      parameter    ( mxspos = 10000,  mxch = 1000 )
c->   REAL           SPOS(MXSPOS,MXCH,3)
c
      real*8         pos(3,mxch), vel(3,mxch)
      real*8         xout(mxch), yout(mxch), zout(mxch)
      real*8         accn(3), delr(3), taccn(3)
c
      integer        others(mxch,mxch)
c
      real*8         atqrp, dt, nu
      integer        nch, nt, savint, trint, plint 
c
      real*8         lsq
      integer        i, ich, it, nchm1, nspos, other, rc
c
      real*8         size
      parameter    ( size = 0.03d0 )
      integer        nside,      colx,     coly,     colz
      parameter    ( nside = 15, colx = 3,  coly = 4,  colz = 6  )
c
      integer        ctlchg / 3 /, ctlpar / 2 /, ctlplt / 4 /,
     *               dstout / 7 /, user / 0 /, posout / 20 /
c
c     Read and evaluate command line.
c
      narg = iargc()
      call getarg(0,pgname)
      if( narg .lt. 1 ) then
999      WRITE(0,*) 'usage: ',PGNAME(1:INDLNB(PGNAME)),
     *              ' <# of charges>'
         stop
      else       
         call getarg(1,argi)
         nch = -1
         read(argi,*) nch
         if( nch .eq. -1  .or.  nch .gt. mxch ) go to 999
      end if
c
c     Open various files.
c
      CALL MYOPEN(CTLPAR,'chgsph.control','OLD',0,0,0,RC)
      if( rc .ne. 0 ) then
         write(0,*) 
     *   '<<< CHGSPH:: Error opening chgsph.control. >>>'
         stop
      else 
c->      WRITE(0,*) 
c->  *   '<<< CHGSPH:: Opened chgsph.control. >>>'
      end if
c->   CALL MYOPEN(CTLCHG,'chgiv.control','OLD',0,0,0,RC)
c->   IF( RC .NE. 0 ) THEN
c->      WRITE(0,*) 
c->  *   '<<< CHGSPH:: Error opening chgiv.control. >>>'
c->      STOP
c->   ELSE 
c->      WRITE(0,*) 
c->  *   '<<< CHGSPH:: Opened chgiv.control. >>>'
c->   END IF
C*    CALL MYOPEN(CTLPLT,'PLOT   CONTROL','OLD',0,0,0,RC)
      CALL MYOPEN(DSTOUT,'chgsph.output',' ',0,0,0,RC)
      if( rc .ne. 0 ) then
         write(0,*) 
     *   '<<< CHGSPH:: Error opening chgsph.output. >>>'
         stop
      else 
c->      WRITE(0,*) 
c->  *   '<<< CHGSPH:: Opened chgsph.output. >>>'
      end if
      CALL MYOPEN(POSOUT,'3d.pos',' ',0,0,0,RC)
      if( rc .ne. 0 ) then
         write(0,*) 
     *   '<<< CHGSPH:: Error opening 3d.pos. >>>'
         stop
      else 
c->      WRITE(0,*) 
c->  *   '<<< CHGSPH:: Opened 3d.pos. >>>'
      end if
c
      read(ctlpar,1010) nt
      read(ctlpar,1000) dt, nu, atqrp
      read(ctlpar,1010) savint, trint, plint 
1000  format(/8f10.3)
1010  format(/8i10)
c
c     Generate random positions velocities.
c    
      do 10 ich = 1 , nch
         pos(1,ich) = 2.0d0 * drand48() - 1.0d0
         pos(2,ich) = 2.0d0 * drand48() - 1.0d0
         pos(3,ich) = 2.0d0 * drand48() - 1.0d0
         vel(1,ich) = 0.0d0
         vel(2,ich) = 0.0d0
         vel(3,ich) = 0.0d0
         call dvsm(pos(1,ich),
     *            1.0d0/dsqrt(dvvdot(pos(1,ich),pos(1,ich),3)),
     *            pos(1,ich),3)
         call dvvprj(vel(1,ich),pos(1,ich),vel(1,ich),3)
 10   continue
c
      nchm1 = nch - 1
      do 30 ich = 1 , nch
         other = 1
         do 20 i = 1 , nch
            if( i .eq. ich ) go to 20
                others(ich,other) = i
                other = other + 1
 20      continue
 30   continue
c
      nspos = 0
c
      do 110 it = 1 , nt
         do 60 ich = 1 , nch
            call dvls(taccn,0.0d0,3)
            do 40 i = 1 , nchm1
               call dvvs(pos(1,ich),pos(1,others(ich,i)),delr,3)
               lsq = dvvdot(delr,delr,3)
               call dvsm(delr,atqrp*1.0d0/(lsq * dsqrt(lsq)),accn,3)
               call dvva(taccn,accn,taccn,3)
 40         continue
            call dvvprj(taccn,pos(1,ich),taccn,3)
            do 50 i = 1 , 3
               taccn(i) = taccn(i) - nu * vel(i,ich)
 50         continue
            call dvsvsa(vel(1,ich),taccn,vel(1,ich),1.0d0,dt,3)
            call dvvprj(vel(1,ich),pos(1,ich),vel(1,ich),3)
            call dvsvsa(pos(1,ich),vel(1,ich),pos(1,ich),1.0d0,dt,3)
            call dvsm(pos(1,ich),
     *                1.0d0/dsqrt(dvvdot(pos(1,ich),pos(1,ich),3)),
     *                pos(1,ich),3)
 60      continue
c->         IF( MYMOD(IT,SAVINT) .NE. 0 ) GO TO 90
c->             NSPOS = NSPOS + 1
c->             DO 80 ICH = 1 , NCH
c->                DO 70 I = 1 , 3
c->                   SPOS(NSPOS,ICH,I) = POS(I,ICH)
c-> 70             CONTINUE
c-> 80          CONTINUE
c-> 90      CONTINUE
         if( mymod(it,trint) .ne. 0 ) go to 100
             write(0,2000) it
2000         FORMAT(' Time step: ',I5)
100      continue
         if( mymod(it,plint) .eq. 0 ) then
c           CALL PLOTIT(POS,NCH,SIZE,NSIDE,COLX,COLY,COLZ)
c           CALL SPHPLOT(POS,NCH,0.05d0)
            do ich = 1 , nch
               xout(ich) = pos(1,ich)
               yout(ich) = pos(2,ich)
               zout(ich) = pos(3,ich)
            end do
            call out_pp3d(it*dt,xout,yout,zout,nch)
         end if
c
110   continue
c   
c     Dump positions.
c
      write(posout,*) nch
      do 2500 ich = 1 , nch
         write(posout,*) pos(1,ich), pos(2,ich), pos(3,ich)
2500  continue
      close(posout)

      if( plint .ne. 0 ) then
c        CALL SPHPLOT(POS,-NCH,0.05d0)
      end if
c
      write(dstout,3000)
3000  FORMAT(' Evolution complete.')
      do 130 ich = 1 , nch
         write(dstout,3010)
3010     format(/)
         do 120 i = 1 , nchm1
            call dvvs(pos(1,ich),pos(1,others(ich,i)),delr,3)
            lsq = dsqrt(dvvdot(delr,delr,3))
            write(dstout,3020) ich, others(ich,i), lsq
3020        FORMAT(' distance(',I2,',',I2,'): ',F12.6)
120      continue
130   continue
c
      stop
c
      end
c
c---------------------------------------------------------------------
c
c     Front end to MOD which prevents division by 0 when second
c     argument is 0, instead returning a non-zero value to
c     expediate I/O decision making etc.
c
c---------------------------------------------------------------------
c
      integer function mymod(i1,i2)
c
         integer     mod
         integer     i1, i2
c
         if( i2 .eq. 0 ) go to 10
             mymod = mod(i1,i2)
         go to 20
 10          mymod = 1
 20      return
c
      end

c---------------------------------------------------------------------
c     Outputs time, charge positions in format suitable for piping 
c     into 'pp3d'.
c---------------------------------------------------------------------
      subroutine out_pp3d(t,x,y,z,n)
         implicit         none

         integer          n
         real*8           t
         real*8           x(n), y(n), z(n)

         real*8           def_size
         parameter      ( def_size = 1.0d0 )

         integer          i,       nfirst
         logical          first
         data             first / .true. /

         save

c---------------------------------------------------------------------
c        Write header info (particle radii)
c---------------------------------------------------------------------
         if( first ) then
            first  = .false.
            nfirst = n
            write(*,*) n
            do i = 1 , nfirst
               write(*,*) def_size
            end do
         end if

         write(*,*) t
         do i = 1 , n
            write(*,*) x(i), y(i), z(i)
         end do

         return

      end 
