function gca_get_particle_dt(partp, end_time) result(dt_p)
use mod_odeint
use mod_global_parameters
type(particle_ptr), intent(inout) :: partp
double precision, intent(in) :: end_time
double precision :: dt_p
double precision :: tout, dt_particles_mype, dt_cfl0, dt_cfl1,&
dt_a
double precision :: dxmin, vp, a, gammap
double precision :: v(ndir), y(ndir+2),ytmp(ndir+2),&
dydt(ndir+2), v0(ndir), v1(ndir), dydt1(ndir+2)
double precision :: ap0, ap1, dt_cfl_ap0, dt_cfl_ap1, dt_cfl_ap
double precision :: dt_euler, dt_tmp
! make these particle cfl conditions more restrictive if you are interpolating out of the grid
double precision :: cfl, uparcfl
double precision, parameter :: uparmin=1.0d-6*const_c
integer :: ipart, iipart, nout, ic1,ic2,ic3,&
igrid_particle, ipe_particle, ipe
logical :: BC_applied
if (const_dt_particles > 0) then
dt_p = const_dt_particles
return
end if
cfl = particles_cfl
uparcfl = particles_cfl
igrid_working = partp%igrid
ipart_working = partp%self%index
dt_tmp = (end_time - partp%self%time)
if(dt_tmp .le. 0.0d0) dt_tmp = smalldouble
! make sure we step only one cell at a time, first check CFL at current location
! then we make an Euler step to the new location and check the new CFL
! we simply take the minimum of the two timesteps.
! added safety factor cfl:
dxmin = min(rnode(rpdx1_,igrid_working),rnode(rpdx2_,igrid_working),&
rnode(rpdx3_,igrid_working),bigdouble)*cfl
! initial solution vector:
y(1:ndir) = partp%self%x(1:ndir) ! position of guiding center
y(ndir+1) = partp%self%u(1) ! parallel momentum component (gamma v||)
y(ndir+2) = partp%self%u(2) ! conserved magnetic moment
ytmp=y
!y(ndir+3) = partp%self%u(3) ! Lorentz factor of guiding centre
call derivs_gca(partp%self%time,y,dydt)
v0(1:ndir) = dydt(1:ndir)
ap0 = dydt(ndir+1)
! guiding center velocity:
v(1:ndir) = abs(dydt(1:ndir))
vp = sqrt(sum(v(:)**2))
dt_cfl0 = dxmin / max(vp, smalldouble)
dt_cfl_ap0 = uparcfl * abs(max(abs(y(ndir+1)),uparmin) / max(abs(ap0),&
smalldouble))
!dt_cfl_ap0 = min(dt_cfl_ap0, uparcfl * sqrt(abs(unit_length*dxmin/(ap0+smalldouble))) )
! make an Euler step with the proposed timestep:
! new solution vector:
dt_euler = min(dt_tmp,dt_cfl0,dt_cfl_ap0)
y(1:ndir+2) = y(1:ndir+2) + dt_euler * dydt(1:ndir+2)
partp%self%x(1:ndir) = y(1:ndir) ! position of guiding center
partp%self%u(1) = y(ndir+1) ! parallel momentum component (gamma v||)
partp%self%u(2) = y(ndir+2) ! conserved magnetic moment
! first check if the particle is outside the physical domain or in the ghost cells
if(.not. particle_in_igrid(ipart_working,igrid_working)) then
y(1:ndir+2) = ytmp(1:ndir+2)
end if
call derivs_gca_rk(partp%self%time+dt_euler,y,dydt)
!call derivs_gca(partp%self%time+dt_euler,y,dydt)
v1(1:ndir) = dydt(1:ndir)
ap1 = dydt(ndir+1)
! guiding center velocity:
v(1:ndir) = abs(dydt(1:ndir))
vp = sqrt(sum(v(:)**2))
dt_cfl1 = dxmin / max(vp, smalldouble)
dt_cfl_ap1 = uparcfl * abs(max(abs(y(ndir+1)),uparmin) / max(abs(ap1),&
smalldouble))
!dt_cfl_ap1 = min(dt_cfl_ap1, uparcfl * sqrt(abs(unit_length*dxmin/(ap1+smalldouble))) )
dt_tmp = min(dt_euler, dt_cfl1, dt_cfl_ap1)
partp%self%x(1:ndir) = ytmp(1:ndir) ! position of guiding center
partp%self%u(1) = ytmp(ndir+1) !parallel momentum component (gamma v||)
partp%self%u(2) = ytmp(ndir+2) ! conserved magnetic moment
!dt_tmp = min(dt_cfl1, dt_cfl_ap1)
! time step due to parallel acceleration:
! The standard thing, dt=sqrt(dx/a) where we compute a from d(gamma v||)/dt and d(gamma)/dt
! dt_ap = sqrt(abs(dxmin*unit_length*y(ndir+3)/( dydt(ndir+1) - y(ndir+1)/y(ndir+3)*dydt(ndir+3) ) ) )
! vp = sqrt(sum(v(1:ndir)**))
! gammap = sqrt(1.0d0/(1.0d0-(vp/const_c)**2))
! ap = const_c**2/vp*gammap**(-3)*dydt(ndir+3)
! dt_ap = sqrt(dxmin*unit_length/ap)
!dt_a = bigdouble
!if (dt_euler .gt. smalldouble) then
! a = sqrt(sum((v1(1:ndir)-v0(1:ndir))**2))/dt_euler
! dt_a = min(sqrt(dxmin/a),bigdouble)
!end if
!dt_p = min(dt_tmp , dt_a)
dt_p = dt_tmp
! Make sure we do not advance beyond end_time
call limit_dt_endtime(end_time - partp%self%time, dt_p)
end function gca_get_particle_dt