---
title: "A Space Charge Calculation"
subtitle: "Example R Notebook"
author: "M J Syphers"
date: "12/23/2016"
output: html_notebook
---


We wish to track a test particle in the presence of the space charge force due to its neighbors (assuming a Gaussian distribution), and determine the *tune* of the particle by performing a Fast Fourier Transform of the tracking data.   By tracking particles with a variety of initial amplitudes, we can see how the tune of the particles varies with amplitude.  We then wish to compare this result to that of the analytical formula found in the textbook by [Edwards and Syphers](http://onlinelibrary.wiley.com/book/10.1002/9783527617272).  We will use the **R** programming environment for the computations.

### The Set-Up

We begin by setting up the problem, defining parameters for the computation.  We'll let $N$ be the number of steps to use for each particle in the calculation.  

Let us assume a space charge tune shift of value $\Delta\nu$.  These parameters are defined here:

```{r setup1}
N = 2^14  # larger number gives better numerical accuracy in the FFT
Dnu = 0.01  # space charge tune shift
```

We next set up arrays for keeping track of a particle's transverse coordinate, $u = x/\sigma$, its associated momentum $du$, and arrays to use in the FFT analysis.

```{r setup2 }
u    <- array(0,N)
du   <- array(0,N)
u0   <- array(0,10)
nu_u <- array(0,10)
```

The unperturbed tune is taken to have value $\nu_0 = 1/(2\pi)$ and hence the phase advance per segment will be 1 radian.  The chunk of code below sets up coefficients used in the tracking.

```{r setup3 }
nu_0 = 1/2/pi
ss = sin(1)
cc = cos(1)
kk = 4*Dnu*2*pi
Amax = 15
```

### Tracking with the Space Charge Force

The space charge "kick" to be coded into the simulation is given by
$$
 \Delta du = 8\pi\Delta\nu\frac{1-e^{-u^2/2}}{u}
$$
which we divide up into two "half-kicks" at the beginning and end of each segment, for symmetry.  The code below loops over $Amax$ = `r Amax` amplitudes and performs an FFT for each amplitude, finding the maximum frequency and storing the result in the array *nu_u*.

```{r tracking }
for(k in 1:Amax){			# Loop over various initial amplitudes
	u0[k] <-  0.2*k
	u[1]  <- u0[k]
	du[1] <- 0 
	for(i in 2:N){
		unxt  <-   u[i-1] ;
		dunxt <-  du[i-1] + kk*(1-exp(-unxt^2/2))/unxt/2 ;		
		u[i]  <-  cc*unxt + ss*dunxt ;
		du[i] <- -ss*unxt + cc*dunxt + kk*(1-exp(-u[i]^2/2))/u[i]/2 ;
	}
	# Fast Fourier analyze the data to extract a "tune"
	#   for Dnu = 0, nu = 1/2/pi
	zz_u <- fft(u)
	zz_u <- abs(zz_u)
	f_u <- c(1:length(u))
	f_u <- f_u/length(u)
	f_ufft <- f_u[ max.col( t( cbind( f_u, zz_u) ) )[2] ]
	nu_u[k] <- min(f_ufft, (1-f_ufft))
}
```

### Plotting the Results

We take a quick look at the data, and see that it has the expected general shape:

```{r quickPlot}
plot(nu_u)
```

To compare with the analytical result in the text, which has the form

$$
\nu-\nu_0 = - \Delta\nu \frac{1-e^{-(u_0/2)^2}I_0((u_0/2)^2)}{(u_0/2)^2},
$$

we need to define a new function in **R**:

```{r u0fcn}
del <- function(u_0){
	-Dnu*( 1- exp(-(u_0/2)^2)*besselI((u_0/2)^2,0))/(u_0/2)^2 }
```

The results of the tracking and the function above can be combined onto the same plot:

```{r plotu0}
plot(u0, nu_u-nu_0)
curve(del(x), add=TRUE)
text( 0.7*max(u0), -Dnu*0.9, 
   expression(-Delta * nu * ~~ frac( 
   (1-e^(-u[0]/2)^2)*~~I[0]((u[0]/2)^2), (u[0]/2)^2 ) ) 
)
```

### Discusssion

From the above, we not only see that the analytical result agrees well with the tracking results, but we also see how straightforward it can be to combine standard LaTeX language and chunks of **R** 
code into a single document, thus creating a much improved method for documenting calculations and results.