c====================================
      subroutine setup1
c====================================
      implicit real*8 (a-h,o-z)

      include 'param.inc'

      common / main1 / x(nox), y(nox), z(nox), is(nox), ic(nox), id
      common / main2 / dif1(nox), dif2(nox)

      common / analyze1 / xdis(-nrg:nrg), xps(3,10,-nrg:nrg), nxpnt
      common / analyze2 / itime, ntime, dt, t0, ntrack
      common / analyze3 / prb(3,nit)
      common / analyze4 / reacd1(3), reacp1(3), ipflg(3), iflg
      common / analyze5 / rdis(0:nrg), rps(3,10,0:nrg), irpoint
      common / analyze6 / reac, scav, escp, DNAe, DNAo, DNAc 


      is(1) = 1             ! 1 = eaq
      is(2) = 2             ! 2 = OH
      is(3) = 3             ! 3 = H3O+

      ic(1) = -1
      ic(2) =  0
      ic(3) =  1

      dif1(1)  =  4.9d-9    ! (m2/s) 
      dif1(2)  =  2.2d-9    ! (m2/s) 
      dif1(3)  =  9.5d-9    ! (m2/s) 


      do i = 1, 3
        dif2(i)  =  dsqrt( 4.d0 * dif1(i) * dt )
      enddo


        xmesh = 0.1d0 * xnm
        nxpnt = 250
      do 10 i = -nxpnt, nxpnt
        xdis(i) = dble(i) * xmesh
      do 10 j = 1, 10
      do 10 k = 1, 3
        xps(k,j,i) = 0.d0
   10 continue


      rmesh  = 1.d-1 * xnm   ! (nm)
      irpoint = 10000
      do 20 i = 0, irpoint 
        rdis(i) = dble(i) * rmesh
      do 20 j = 1, 10
      do 20 k = 1, 3
        rps(k,j,i) = 0.d0
   20 continue



      reacd1(2) = 0.72d-9        ! (m)
      reacd1(3) = 0.75d-9        ! (m)

      reacp1(2) = 4.9d-1
      reacp1(3) = 3.8d-2


      do it = 0, itime
      do i = 1, 3
        prb(i,it) = 1.d0
      enddo
      enddo


      reac = 0.d0
      scav = 0.d0
      escp = 0.d0
      DNAe = 0.d0
      DNAo = 0.d0
      DNAc = 0.d0

      end      


c====================================
      subroutine setup2
c====================================
      implicit real*8 (a-h,o-z)

      include 'param.inc'

      common / main1 / x(nox), y(nox), z(nox), is(nox), ic(nox), id

      common / analyze4 / reacd1(3), reacp1(3), ipflg(3), iflg
      common / analyze7 / iflge, iflgo, iflgc 

      common / randam / idum

      dimension rx(0:5000), enumg(0:5000), enume(0:5000), enum(0:5000)

      real ran2


*(selection of distribution type)
      if( id .eq. 1 ) goto 11   !  delta
      if( id .eq. 2 ) goto 12   !  gaussian
      if( id .eq. 3 ) goto 13   !  exponential
      if( id .eq. 4 ) goto 14   !  dmcc_phys


*(delta)
   11 continue
        xint = 0.97d0
        goto 15


*(gaussian)
   12 continue

        do i = 0, 5000
          rx(i)     =  0.d0
          enumg(i)  =  0.d0
        enddo

        dr0 = 1.d-2
        n   = 0

        r0    =  8.00d0
        sgm   =  r0 / dsqrt(8.d0/pi)
        fct0  =  dsqrt(8.d0 * pi**3 * sgm**6)

      do r = dr0, 20.d0, dr0
        rx(n+1)  =  r
        fac0     =  r**2 / 2.d0 / sgm**2
        fg       =  4.d0 * pi * r**2 * dexp( -fac0 ) / fct0 * dr0
        enumg(n+1) =  enumg(n) + fg

c        write(*,100) rx(n+1), fg, enumg(n+1)

        n = n + 1
      enddo
c      stop


      rnd = dble( ran2(idum) )
      do i = 0, 2000   ! 525
        if( rnd .ge. enumg(i) .and. rnd .lt. enumg(i+1) ) then
          xint = rx(i) + dr0 / 2.d0
        endif
      enddo
      goto 15


*(expornential) 
   13 continue

        do i = 0, 5000
          rx(i)     =  0.d0
          enume(i)  =  0.d0
        enddo

        dr0  =  1.d-2
        n   =  0

        r1    =  4.50d0
        b     =  r1 / 3.d0
        fct1  =  8.d0 * pi * b**3

      do r = dr0, 15.d0, dr0  
        rx(n+1)  =  r
        fac1     =  r / b
        fe       =  4.d0 * pi * r**2 * dexp( -fac1 ) / fct1 * dr0   
        enume(n+1) =  enume(n) + fe 

c      write(*,100) rx(n+1), fe, enume(n+1)      

        n = n + 1
      enddo 

c      stop


      rnd = dble( ran2(idum) )
      do i = 0, 1500  ! 525
        if( rnd .ge. enume(i) .and. rnd .lt. enume(i+1) ) then
          xint = rx(i) + dr0 * dble( ran2(idum) ) 
        endif
      enddo
      goto 15


*(dmcc_phys)
   14 continue


         iradio = 0

      if( iradio .eq. 1 ) then

           dr0 = 0.1d0

         do i = 0, 1001
           rx(i)    =  0.d0
           enum(i)  =  0.d0
         enddo

        open(1,file='rdist.out')
         do i = 0, 1000
           read(1,*) rx(i), res1, res2, res3, res4, res5, 
     .                      res6, res7, res8, res9, res10
           if( i .eq. 0 ) then
             enum(i) = res6
           else
             enum(i) = enum(i-1) + res6
           endif
c        write(*,100) rx(i), enum(i), 1.d0-enum(i)
         enddo
        close(1)
c        stop

        rnd = dble( ran2(idum) )
        do i = 0, 1000
          if( rnd .ge. enum(i) .and. rnd .lt. enum(i+1) ) then
            xint = rx(i) + dr0 * dble( ran2(idum) ) 
          endif
        enddo


      else

           dr0 = 0.1d0

         do i = 0, 251
           rx(i)    =  0.d0
           enum(i)  =  0.d0
         enddo

        open(1,file='Sout.dat')
         do i = 1, 251
           read(1,*) rx(i), res1, res2, res3, res4, res5, res6
           enum(i) = enum(i-1) + res6
         enddo
        close(1)


        rnd = dble( ran2(idum) )
        do i = 1, 251
          if( rnd .ge. enum(i) .and. rnd .lt. enum(i+1) ) then
            xint = rx(i) + dr0 * dble( ran2(idum) )
          endif
        enddo

      endif

  100 format(5(1pe15.7))
   15 continue


        the = dble( ran2(idum) ) * pi
        phi = dble( ran2(idum) ) * pi * 2.d0
        r   = xint
        sft = 6.00d0
        
        x(1) = ( r * dsin( the ) * dcos( phi ) + sft ) * xnm  ! (m)
        y(1) = ( r * dsin( the ) * dsin( phi )       ) * xnm  ! (m)
        z(1) = ( r * dcos( the )                     ) * xnm  ! (m)

c        x(1) = xint * xnm  ! (m)
c        y(1) = 0.d0
c        z(1) = 0.d0


      do i = 2, 3
        x(i) = sft * xnm   ! (m)
        y(i) = 0.d0
        z(i) = 0.d0
      enddo


        ipflg(2) = 0
        ipflg(3) = 0
        iflg     = 0
        iflge    = 0
        iflgo    = 0
        iflgc    = 0
 
      end


c====================================
      subroutine analyze( it ) 
c====================================
      implicit real*8 (a-h,o-z)

      include 'param.inc'

      common / main1 / x(nox), y(nox), z(nox), is(nox), ic(nox), id
      common / main3 / dx(3,3), dy(3,3), dz(3,3), dr(3,3)

      common / analyze1 / xdis(-nrg:nrg), xps(3,10,-nrg:nrg), nxpnt
      common / analyze2 / itime, ntime, dt, t0, ntrack
      common / analyze5 / rdis(0:nrg), rps(3,10,0:nrg), irpoint

      dimension out(10)


      t = t0 + dble( it - 1 ) * dt


      out( 1) = 5.d-12 
      out( 2) = 3.d-11 
      out( 3) = 5.d-11 
      out( 4) = 1.d-10 
      out( 5) = 3.d-10 
      out( 6) = 5.d-10 
      out( 7) = 1.d-9 
      out( 8) = 3.d-9 
      out( 9) = 5.d-9 
      out(10) = 1.d-8 


c      do itd = 1, 7
c        if( t .eq. out(itd) ) then
c          do ii = 1, 3
c            do i = -nxpnt, nxpnt
c              if( x(ii).ge.xdis(i) .and. x(ii).lt.xdis(i+1) ) then
c                xps(ii,itd,i)  =  xps(ii,itd,i) + 1.d0 / dble(ntrack)
c              endif
c            enddo
c          enddo
c        endif
c      enddo 


      do itd = 1, 10
        if( t .eq. out(itd) ) then
          do ii = 2, 3
            do i = 0, irpoint
              if( dr(1,ii).ge.rdis(i) .and. dr(1,ii).lt.rdis(i+1) ) then
                 rps(ii,itd,i)  =  rps(ii,itd,i) + 1.d0 / dble(ntrack)
              endif
            enddo
          enddo
        endif
      enddo


c      xe =  x(1)
c      ye =  y(1)
c      ze =  z(1)
c      rr =  dsqrt( xe**2 + ye**2 + ze**2 )


c      do itd = 1, 4 
c        if( t .eq. out(itd) ) then
c          do ii = 1, 1
c            do i = 0, irpoint
c              if( rr.ge.rdis(i) .and. rr.lt.rdis(i+1) ) then
c                 rps(ii,itd,i)  =  rps(ii,itd,i) + 1.d0 / dble(ntrack) 
c              endif
c            enddo
c          enddo
c        endif 
c      enddo

      end


c===========================================
      subroutine output(fname1, fname2)
c===========================================
      implicit real*8 (a-h,o-z)

      include 'param.inc'

      common / analyze1 / xdis(-nrg:nrg), xps(3,10,-nrg:nrg), nxpnt
      common / analyze2 / itime, ntime, dt, t0, ntrack
      common / analyze3 / prb(3,nit)
      common / analyze5 / rdis(0:nrg), rps(3,10,0:nrg), irpoint
      common / analyze6 / reac, scav, escp, DNAe, DNAo, DNAc

      character*12 fname1, fname2


* (scale)
       ps  =  1.d-12            ! (pico-second) 


      do i = -nxpnt, nxpnt
c        write(*,100) xdis(i)/xnm, (xps(1,itd,i),itd=1,3)
      enddo


      open(1,file= fname1 )
        do i = 0, irpoint 
          write(1,100) rdis(i)/xnm, ( rps(2,itd,i), itd= 1, 10 )
        enddo
      close(1)


      open(1,file= fname2 )
        do i = 0, irpoint
          write(1,100) rdis(i)/xnm, ( rps(3,itd,i), itd= 1, 10 )
        enddo
      close(1)


      do it = 1, itime

        t  =  ( t0 + dble(it-1) * dt ) 

          
        if( it .eq. 1 .or. mod(it,ntime) .eq. 0 ) then
          write(*,100) t/ps, ( prb(i,it) / prb(i,1), i = 1, 3 ) 
        endif

c        if( mod(it,ntime) .eq. 0 ) then
c          write(*,100) t*1.d+3, prb(1,it), prb(1,it)*6.22d0
c        endif

      enddo
 

      write(*,100) escp, reac, scav, DNAe, DNAo, DNAc


  100 format(15(1pe15.7))
      end


