*Originally Posted by ***b49P23TIvg**

itemindex=numpy.where(array==item)

If i and j are determined by

numpy.where(array==item)

I do not understand this part.
*Originally Posted by ***b49P23TIvg**

polyfunc = polyfit(s.t[i:j],s.firstharmonic[i:j],1)

That would mean that I would have to always explicitly know what i and j are. The Two Stream instability model isn't the same everytime because it uses generators. Its like my previous Landau Damping but in reverse.
Code:

#!/usr/bin/env python
#
# Electrostatic PIC code in a 1D cyclic domain
from numpy import *
from scipy import *
import matplotlib.pyplot as plt # Matplotlib plotting library
try:
import matplotlib.gridspec as gridspec # For plot layout grid
got_gridspec = True
except:
got_gridspec = False
# Need an FFT routine, either from SciPy or NumPy
try:
from scipy.fftpack import fft, ifft
except:
# No SciPy FFT routine. Import NumPy routine instead
from fft import fft, ifft
def rk4step(f, y0, dt, args=()):
""" Takes a single step using RK4 method """
k1 = f(y0, *args)
k2 = f(y0 + 0.5*dt*k1, *args)
k3 = f(y0 + 0.5*dt*k2, *args)
k4 = f(y0 + dt*k3, *args)
return y0 + (k1 + 2.*k2 + 2.*k3 + k4)*dt / 6.
def calc_density(position, ncells, L):
""" Calculate charge density given particle positions
Input
position - Array of positions, one for each particle
assumed to be between 0 and L
ncells - Number of cells
L - Length of the domain
Output
density - contains 1 if evenly distributed
"""
# This is a crude method and could be made more efficient
density = zeros([ncells])
nparticles = len(position)
dx = L / ncells # Uniform cell spacing
for p in position / dx: # Loop over all the particles, converting position into a cell number
plower = int(p) # Cell to the left (rounding down)
offset = p - plower # Offset from the left
density[plower] += 1. - offset
density[(plower + 1) % ncells] += offset
# nparticles now distributed amongst ncells
density *= float(ncells) / float(nparticles) # Make average density equal to 1
return density
def periodic_interp(y, x):
"""
Linear interpolation of a periodic array y at index x
Input
y - Array of values to be interpolated
x - Index where result required. Can be an array of values
Output
y[x] with non-integer x
"""
ny = len(y)
if len(x) > 1:
y = array(y) # Make sure it's a NumPy array for array indexing
xl = floor(x).astype(int) # Left index
dx = x - xl
xl = ((xl % ny) + ny) % ny # Ensures between 0 and ny-1 inclusive
return y[xl]*(1. - dx) + y[(xl+1)%ny]*dx
def fft_integrate(y):
""" Integrate a periodic function using FFTs
"""
n = len(y) # Get the length of y
f = fft(y) # Take FFT
# Result is in standard layout with positive frequencies first then negative
# n even: [ f(0), f(1), ... f(n/2), f(1-n/2) ... f(-1) ] Only want the first and the last
# n odd: [ f(0), f(1), ... f((n-1)/2), f(-(n-1)/2) ... f(-1) ]
"""Just interested in the first harmonic for the fastest growing mode"""
if n % 2 == 0: # If an even number of points (if not then n is odd)
k = concatenate( (arange(0, n/2+1), arange(1-n/2, 0)) )
else:
k = concatenate( (arange(0, (n-1)/2+1), arange( -(n-1)/2, 0)) )
k = 2.*pi*k/n
# Modify frequencies by dividing by ik
f[1] /= (1j * k[1])
f[n-1] /= (1j * k[n-1]) #Have cut out alot of harmonics here
f[2:n] = 0.0
f[0] = 0. # Set the arbitrary zero-frequency term to zero
return ifft(f).real # Reverse Fourier Transform
def pic(f, ncells, L):
""" f contains the position and velocity of all particles
"""
nparticles = len(f)/2 # Two values for each particle
pos = f[0:nparticles] # Position of each particle
vel = f[nparticles:] # Velocity of each particle
dx = L / float(ncells) # Cell spacing
# Ensure that pos is between 0 and L
pos = ((pos % L) + L) % L
# Calculate number density, normalised so 1 when uniform
density = calc_density(pos, ncells, L)
# Subtract ion density to get total charge density
rho = density - 1.
# Calculate electric field
E = -fft_integrate(rho)*dx
# Interpolate E field at particle locations
accel = -periodic_interp(E, pos/dx)
# Put back into a single array
return concatenate( (vel, accel) )
####################################################################
def run(pos, vel, L, ncells=None, out=[], output_times=linspace(0,20,100), cfl=0.5):
if ncells == None:
ncells = int(sqrt(len(pos))) # A sensible default
dx = L / float(ncells)
f = concatenate( (pos, vel) ) # Starting state
nparticles = len(pos)
time = 0.0
for tnext in output_times:
# Advance to tnext
stepping = True
while stepping:
# Maximum distance a particle can move is one cell
dt = cfl * dx / max(abs(vel))
if time + dt >= tnext:
# Next time will hit or exceed required output time
stepping = False
dt = tnext - time
#print "Time: ", time, dt
f = rk4step(pic, f, dt, args=(ncells, L))
time += dt
# Extract position and velocities
pos = ((f[0:nparticles] % L) + L) % L
vel = f[nparticles:]
# Send to output functions
for func in out:
func(pos, vel, ncells, L, time)
return pos, vel
####################################################################
#
# Output functions and classes
#
class Plot:
"""
Displays three plots: phase space, charge density, and velocity distribution
"""
def __init__(self, pos, vel, ncells, L):
d = calc_density(pos, ncells, L)
vhist, bins = histogram(vel, int(sqrt(len(vel))))
vbins = 0.5*(bins[1:]+bins[:-1])
# Plot initial positions
if got_gridspec:
self.fig = plt.figure()
self.gs = gridspec.GridSpec(4, 4)
ax = self.fig.add_subplot(self.gs[0:3,0:3])
self.phase_plot = ax.plot(pos, vel, '.')[0]
ax.set_title("Phase space")
ax = self.fig.add_subplot(self.gs[3,0:3])
self.density_plot = ax.plot(linspace(0, L, ncells), d)[0]
ax = self.fig.add_subplot(self.gs[0:3,3])
self.vel_plot = ax.plot(vhist, vbins)[0]
else:
self.fig = plt.figure()
self.phase_plot = plt.plot(pos, vel, '.')[0]
self.fig = plt.figure()
self.density_plot = plt.plot(linspace(0, L, ncells), d)[0]
self.fig = plt.figure()
self.vel_plot = plt.plot(vhist, vbins)[0]
plt.ion()
plt.show()
def __call__(self, pos, vel, ncells, L, t):
d = calc_density(pos, ncells, L)
vhist, bins = histogram(vel, int(sqrt(len(vel))))
vbins = 0.5*(bins[1:]+bins[:-1])
self.phase_plot.set_data(pos, vel) # Update the plot
self.density_plot.set_data(linspace(0, L, ncells), d)
self.vel_plot.set_data(vhist, vbins)
plt.draw()
class Summary:
def __init__(self):
self.t = []
self.firstharmonic = []
def __call__(self, pos, vel, ncells, L, t):
# Calculate the charge density
d = calc_density(pos, ncells, L)
# Amplitude of the first harmonic
fh = 2.*abs(fft(d)[1]) / float(ncells)
print "Time:", t, "First:", fh
self.t.append(t)
self.firstharmonic.append(fh)
####################################################################
#
# Functions to create the initial conditions
#
def landau(npart, L, alpha=0.2):
"""
Creates the initial conditions for Landau damping
"""
# Start with a uniform distribution of positions
pos = random.uniform(0., L, npart)
pos0 = pos.copy()
k = 2.*pi / L
for i in range(10): # Adjust distribution using Newton iterations
pos -= ( pos + alpha*sin(k*pos)/k - pos0 ) / ( 1. + alpha*cos(k*pos) )
# Normal velocity distribution
vel = random.normal(0.0, 1.0, npart)
return pos, vel
def twostream(npart, L, vbeam=2):
# Start with a uniform distribution of positions
pos = random.uniform(0., L, npart)
# Normal velocity distribution
vel = random.normal(0.0, 1.0, npart)
np2 = int(npart / 2)
vel[:np2] += vbeam # Half the particles moving one way
vel[np2:] -= vbeam # and half the other
return pos,vel
####################################################################
if __name__ == "__main__":
#vlist = [1., 5., 10.]
# 2-stream instability
#for v in vlist:
L = 100
ncells = 20
pos, vel = twostream(10000, L, 3.)
# Create some output classes
#p = Plot(pos, vel, ncells, L) # This displays an animated figure
s = Summary() # Calculates, stores and prints summary info
# Run the simulation
pos, vel = run(pos, vel, L, ncells,
out=[s],
#out=[p, s], # These are called each output
output_times=linspace(0.,40,100)) # The times to output
x = s.t #time
data = s.firstharmonic #harmonics y values
c = (diff(sign(diff(data))) < 0).nonzero()[0] + 1 # local max
timemax = []
maxh = []
rmaxh = []
for i in c:
timemax.append(x[i])
maxh.append(data[i])
#print maxh
for i in reversed(maxh):
print i
rmaxh.append(i)
#print rmaxh
def crap(rmaxh):
for j in range(len(rmaxh)-1):
if not rmaxh[j] > rmaxh[j+1]:
sliceindex = j
slicevalue = rmaxh[j]
return sliceindex
sliceindex = crap(rmaxh)
def twostream(timemax, rmaxh):
timestream = timemax[: sliceindex]
timenoisy = timemax[sliceindex :]
rmaxhstream = rmaxh[: sliceindex]
rmaxhnoisy= rmaxh[sliceindex:]
print timestream
print timenoisy
#coefs = polyfit(timestream, maxhstream,1)
#straightfit = polyval(polyfunc,timestream)
#return timestream, straightfit
#timestream, straightfit, = twostream(timemax, log(rmaxh))
# Summary stores an array of the first-harmonic amplitude
# Make a semilog plot to see exponential damping
plt.figure()
plt.plot(s.t, s.firstharmonic)
plt.plot(timemax, maxh, 'or')
plt.xlabel("Time [Normalised]")
plt.ylabel("First harmonic amplitude [Normalised]")
#plt.yscale('log')
#plt.legend()
plt.ioff() # This so that the windows stay open
plt.show()

I have not included our revised technique for the density function because this is a fresh code..
Tweet This+ 1 thisPost To Linkedin