       integer   maxshells
       real*8    max_p
c
       parameter (maxshells = 200, max_p = 10.0D0) 
c
       real*8      R_shell_x(maxshells), f_acc, temp   
       real*8      num_ptcls_tot, shell_vol, diff_shell_vol, tot_vol
       real*8      num_ptcls, vel, energy, tot_ptcls
       real*8      vo, mp, injectrate
       real*8      del_p, little_q, p 
       real*8      shock_speed, Machno, comp_rat        
       real*8      prev_shock_posn, prev_shock_time 
       real*8      p_inj, p_cut, distr_array(maxshells, 250)     
       real*8      esc_shell(in, 250)                    
       real*8      esc_dens_up(250), esc_dens_down(250),
     *             delta, prev_acc_time, prev_expand_time(maxshells)
       real*8      kappa(250, 1000)
c
       integer   shock_begin, shock_end, shell_posn 
       integer   i , j , k, R_shell_i(maxshells), no_shells
      subroutine add_acc_ptcls
      integer   max_esc_shell_j, p_int
c
       common /shell/  esc_shell,R_shell_x,R_shell_i,no_shells
       common /shock1/ prev_shock_posn, prev_shock_time
       common /shock2/ shock_begin, shock_end, shock_speed,
     *                 Machno, comp_rat 
       common /ptcls/  p_inj, p_cut, distr_array, esc_dens_up,
     *                 esc_dens_down, prev_acc_time, prev_expand_time
       common /diff_coeff/ kappa   
c                       

      max_esc_shell_j = 0   ! Max array variable of esc shell - still trapped.   
c
       vo = 52483.25D0
       mp = 1.67D-27
c
       del_p = (max_p - 2.0)/250.0
c                                         
       little_q = 3.0*comp_rat/(comp_rat - 1.0)
c       
       delta = 0.01D0
       injectrate = delta*d(shock_end+5,1,5)*(shock_speed -   ! d is the density.
     *              v1(shock_end+5,1,5))
c                        
       p = 2.0  ! Momentum is stored logarithmically in the array.              
c 
       tot_ptcls = 0.0D0                         
c
       Do i = 1, 250  
c                        
         p_int = INT((p - 2.0)/del_p + 0.5) + 1  ! Momentum array index.     
c              
         If ((p .GE. Log(p_inj)) .AND. 
     *            (p .LE. Log(p_cut))) Then
c
           f_acc = little_q*injectrate/4.0D0/pi/p_inj**3.0
	   f_acc = f_acc/(shock_speed-v1(shock_end+5,1,5))
	   f_acc = f_acc*(Exp(p)/p_inj)**(-1.0*little_q)
c     
         Else
c
           f_acc = 0.0 
c
         End If 
c
c         tot_ptcls = tot_ptcls + f_acc*4.0*pi*Exp(p)**2.0D0*
c     *               (Exp(p+del_p/2.0D0)-Exp(p-del_p/2.0D0))	 
c                         
         shell_vol = 4.0*3.14159*
     *               R_shell_x(no_shells+1)**3.0/3.0     ! Shell volume.
         shell_vol = shell_vol - 
     *        4.0*3.14159*R_shell_x(no_shells)**3.0/3.0 
c                        
         vel = Exp(p)*52483.25/Sqrt(1.0 + Exp(p)*Exp(p)*
     *         52483.25*52483.25/(3.0E8*3.0E8))          
c                         
         Energy = Sqrt((1.67E-27*3.0E8**2.0)**2.0 + 
     *            (1.67E-27*Exp(p)*52483.25*3.0E8)**2.0) -
     *            1.67E-27*3.0E8**2.0

         Energy = Energy/1.6E-19 
c  
         temp = kappa(i,1)/(shock_speed-v1(shock_end+2,1,5))/vo
	 temp = temp/1.5D11
c               
         num_ptcls_tot = 0.0
         tot_vol = 0.0
c                         
         diff_shell_vol = shell_vol + 4.0*3.14159*        ! Eff. volume that
     *        R_shell_x(no_shells+1)**2.0*temp            ! conserves ptcl
         diff_shell_vol = diff_shell_vol + 8.0*3.14159*   ! number.
     *        R_shell_x(no_shells+1)*temp**2.0 +
     *        8.0*3.14159*temp**3.0
c                         
         Do j = R_shell_i(no_shells)+1, in
c
           num_ptcls = 0.0
c
           If (j .LE. shock_begin) Then    
c
             If ((j-1) .EQ. R_shell_i(no_shells)) Then
c                        
               tot_vol = tot_vol + 
     *                 4.0*3.14159*(x1a(j)**3.0 - 
     *                 R_shell_x(no_shells)**3.0)/3.0                        
c                
               num_ptcls = f_acc*         
     *                   4.0*3.14159*(x1a(j)**3.0 - 
     *                   R_shell_x(no_shells)**3.0)/3.0  
c                                         
             Else
c
               tot_vol = tot_vol +
     *                 4.0*3.14159*(x1a(j)**3.0 - x1a(j-1)**3.0)/3.0
c
               num_ptcls = f_acc*
     *                 4.0*pi*(x1a(j)**3.0 - x1a(j-1)**3.0)/3.0
c
             End If	                         
c
             tot_ptcls = tot_ptcls + num_ptcls*4.0D0*pi*Exp(p)**2.0*
     *                   (Exp(p+del_p/2.0D0)-Exp(p-del_p/2.0D0))
c
             num_ptcls_tot = num_ptcls_tot + num_ptcls            
c
             distr_array(no_shells, p_int) = 
     *          distr_array(no_shells, p_int) + num_ptcls
c                                                     
           Else If ((x1a(j) .GT. x1a(shock_begin)) .AND.     
     *         (x1a(j) .LE. (x1a(shock_begin) + 4.0*temp))) Then        
c    
c             call esc_distr(num_ptcls, f_acc, diff_shell_vol, 
c     *                temp, x1a(j-1), x1a(j))
c
             num_ptcls = f_acc*DExp(-1.0D0*(x1a(j)-x1a(shock_begin))/
     *                               temp)
             num_ptcls = num_ptcls*4.0D0*pi*
     *                (x1a(j)**3.0D0 - x1a(j-1)**3.0D0)/3.0D0      	     
c                                                
             esc_shell(j, p_int) = esc_shell(j, p_int) +
     *               num_ptcls
c
c             distr_array(no_shells, p_int) =
c     *          distr_array(no_shells,p_int) + num_ptcls
c
             tot_ptcls = tot_ptcls + num_ptcls*4.0D0*pi*Exp(p)**2.0*
     *                   (Exp(p+del_p/2.0D0)-Exp(p-del_p/2.0D0))
c
             num_ptcls_tot = num_ptcls_tot + num_ptcls 
c
           Else If ((x1a(j) .GT. (x1a(shock_begin) + 4.0*temp)) .AND.
     *         (x1a(j - 1) .LT. (x1a(shock_begin) + 4.0*temp))) Then 
c
             num_ptcls = f_acc*DExp(-1.0D0*(x1a(j)-x1a(shock_begin))/
     *                               temp)
             num_ptcls = num_ptcls*4.0D0*pi*
     *                   ((x1a(shock_begin)+4.0*temp)**3.0D0-
     *                   x1a(j-1)**3.0D0)/3.0D0
c
             esc_shell(j, p_int) = esc_shell(j, p_int) +
     *             num_ptcls
c
c             distr_array(no_shells, p_int) =
c     *          distr_array(no_shells, p_int) + num_ptcls
c
             tot_ptcls = tot_ptcls + num_ptcls*4.0D0*pi*Exp(p)**2.0*
     *                   (Exp(p+del_p/2.0D0)-Exp(p-del_p/2.0D0))
c
             num_ptcls_tot = num_ptcls_tot + num_ptcls
c	     
             esc_dens_up(p_int) = esc_dens_up(p_int) +
     *                  f_acc*DExp(-1.0D0*(4.0*temp)/temp)  
c                                                     
           End If   
c                                                     
         End Do   
c
         If ((Exp(p) .GE. p_inj) .and. (Exp(p) .LE. p_cut)) Then
	   Write(38,*) Energy, esc_dens_up(p_int), 
     *                 x1a(shock_begin)+4.0D0*temp, time
         Else
	   Write(38,*) Energy, esc_dens_up(p_int),
     *                 x1a(shock_begin), time
         EndIf
         Write(46,*) Energy, f_acc, esc_dens_up(p_int)
c      
         p = p + del_p          
c                                                                         
       End Do 
c                                   
       return
       end
