particleinterp: Interpolate particle coordinates between snapshots

In obreschkow/simstar: Reading, processing and visualising astrophysical N-body simulation data

Description

Interpolate/extrapolate particle positions and velocities between discrete time steps. An independent polynomial interpolation along each Cartesian coordinate is used by default. If an acceleration field is specified via the optional argument afield, a leapfrog integrator with the custom time step dt is applied. This is particularly useful for interpolation between time steps separated by one or more characteristic orbital times.

Usage

particleinterp(
  x0,
  x1,
  t0,
  t1,
  ti = 0.5,
  v0 = NULL,
  v1 = NULL,
  afield = NULL,
  dt = NULL
)

Arguments

x0	n-by-d array of n particle positions in d dimensions at time t0; or an n-vector for one dimension
x1	n-by-d array of n particle positions in d dimensions at time t1; or an n-vector for one dimension
t0	single number giving the time of the first time step, corresponding to the particle positions x0
t1	single number giving the time of the second time step, corresponding to the particle positions x1
ti	single number giving the time of the interpolated/extrapolated output positions and velocities
v0	optional n-by-d array of n particle velocities in d dimensions at time t0; or an n-vector for one dimension
v1	optional n-by-d array of n particle velocities in d dimensions at time t1; or an n-vector for one dimension
afield	optional acceleration field for leapfrog integration. This must be a function of a n-by-d array representing n position vectors, which returns an n-by-d array with the n acceleration vectors in the same length and time units as x and v.
dt	optional time step used for leapfrog integration; only used if afield given. The default is dt=abs(t1-t0)/20.

Author(s)

Danail Obreschkow

Examples

## Example 1: 1D-harmonic oscillator

# analytical model of continuous orbit in a 1D harmonic oscillator
x = function(t) sin(t) # position
v = function(t) cos(t) # velocity
curve(x,0,pi,200,xlab='Time t',ylab='Position x(t)')

# evaluation at two discrete times
t0 = 0.8; t1 = 2.8
points(c(t0,t1),x(c(t0,t1)),pch=16)

# interpolation/extrapolation
for (ti in seq(0,5,by=0.1)) {
  
  # using only positions (blue)
  p = particleinterp(x(t0),x(t1),t0,t1,ti)
  points(ti,p$x,col='blue')
  
  # using positions and velocities (red)
  p = particleinterp(x(t0),x(t1),t0,t1,ti,v(t0),v(t1))
  points(ti,p$x,col='red')
  
}

## Example 2: particles moving in a spherical gravitational potential

# make acceleration field of an isothermal potential of unit circular velocity
afield = function(x) -x/rowSums(x^2)

# analytic function of circular orbits
x = function(t) cbind(x=r*cos(t/r+phi),y=r*sin(t/r+phi)) # positions
v = function(t) cbind(-sin(t/r+phi),cos(t/r+phi)) # velocities 

# specification of two circular orbits
r = c(1,2) # radii of circular orbits
phi = c(0,2) # azimuths at t=0 in radians

# times of the two time steps
t0 = 0 # initial time
t1 = 6 # final time

# plot initial and final positions as black points
plot(rbind(x(t0),x(t1)),pch=16,asp=1,xlim=c(-3,3),ylim=c(-3,3))

# draw lines for analytical circular orbits
plotrix::draw.arc(0,0,r,t0/r+phi,t1/r+phi)

# plot leapfrog-interpolated points as red crosses
for (ti in seq(t0,t1,length=20)) {
  out = particleinterp(x(t0),x(t1),t0,t1,ti,v(t0),v(t1),afield)
  points(out$x,col='red',pch=4)
}

# plot standard polynomial interpolation as blue circles
for (ti in seq(t0,t1,length=20)) {
  out = particleinterp(x(t0),x(t1),t0,t1,ti,v(t0),v(t1))
  points(out$x,col='blue')
}

obreschkow/simstar documentation built on Jan. 29, 2022, 2:16 p.m.