diff --git a/examples/Feedback/feedback.pro b/examples/Feedback/feedback.pro
deleted file mode 100644
index 02d616fc82f0aeb7011d022d13db9d1d1030e89c..0000000000000000000000000000000000000000
--- a/examples/Feedback/feedback.pro
+++ /dev/null
@@ -1,24 +0,0 @@
-base = 'Feedback'
-inf = 'Feedback_005.hdf5'
-
-blast = [5.650488e-01, 5.004371e-01, 5.010494e-01] ; location of blast
-pos = h5rd(inf,'PartType0/Coordinates')
-vel = h5rd(inf,'PartType0/Velocities')
-rho = h5rd(inf,'PartType0/Density')
-utherm = h5rd(inf,'PartType0/InternalEnergy')
-
-; shift to centre
-for ic=0,2 do pos[ic,*] = pos[ic,*] - blast[ic]
-
-;; distance from centre
-dist = fltarr(n_elements(rho))
-for ic=0,2 do dist = dist + pos[ic,*]^2
-dist = sqrt(dist)
-
-; radial velocity
-vr = fltarr(n_elements(rho))
-for ic=0,2 do vr = vr + pos[ic,*]*vel[ic,*]
-vr = vr / dist
-
-;
-end
diff --git a/examples/Feedback/feedback.yml b/examples/Feedback/feedback.yml
deleted file mode 100644
index de4f7abef1ef538a97a5e38c72b4db5ce2647976..0000000000000000000000000000000000000000
--- a/examples/Feedback/feedback.yml
+++ /dev/null
@@ -1,43 +0,0 @@
-# Define the system of units to use internally.
-InternalUnitSystem:
- UnitMass_in_cgs: 1 # Grams
- UnitLength_in_cgs: 1 # Centimeters
- UnitVelocity_in_cgs: 1 # Centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
-
-# Parameters governing the time integration
-TimeIntegration:
- time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 5e-2 # The end time of the simulation (in internal units).
- dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
- dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
-
-# Parameters governing the snapshots
-Snapshots:
- basename: Feedback # Common part of the name of output files
- time_first: 0. # Time of the first output (in internal units)
- delta_time: 1e-2 # Time difference between consecutive outputs (in internal units)
-
-# Parameters governing the conserved quantities statistics
-Statistics:
- delta_time: 1e-3 # Time between statistics output
-
-# Parameters for the hydrodynamics scheme
-SPH:
- resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
- delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
- CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
-
-# Parameters related to the initial conditions
-InitialConditions:
- file_name: ./Feedback.hdf5 # The file to read
-
-# Parameters for feedback
-
-SN:
- time: 0.001 # time the SN explodes (internal units)
- energy: 1.0 # energy of the explosion (internal units)
- x: 0.5 # x-position of explostion (internal units)
- y: 0.5 # y-position of explostion (internal units)
- z: 0.5 # z-position of explostion (internal units)
diff --git a/examples/Feedback/makeIC.py b/examples/Feedback/makeIC.py
deleted file mode 100644
index bd1081a9c275616038f5fa4e3eb943c36cb4c3eb..0000000000000000000000000000000000000000
--- a/examples/Feedback/makeIC.py
+++ /dev/null
@@ -1,109 +0,0 @@
-###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
- # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- # 2016 Tom Theuns (tom.theuns@durham.ac.uk)
- #
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU Lesser General Public License as published
- # by the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- #
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- # GNU General Public License for more details.
- #
- # You should have received a copy of the GNU Lesser General Public License
- # along with this program. If not, see <http://www.gnu.org/licenses/>.
- #
- ##############################################################################
-
-import h5py
-import sys
-from numpy import *
-
-# Generates a swift IC file containing a cartesian distribution of particles
-# at a constant density and pressure in a cubic box
-
-# Parameters
-periodic= 1 # 1 For periodic box
-boxSize = 1.
-L = int(sys.argv[1]) # Number of particles along one axis
-rho = 1. # Density
-P = 1.e-6 # Pressure
-gamma = 5./3. # Gas adiabatic index
-eta = 1.2349 # 48 ngbs with cubic spline kernel
-fileName = "Feedback.hdf5"
-
-#---------------------------------------------------
-numPart = L**3
-mass = boxSize**3 * rho / numPart
-internalEnergy = P / ((gamma - 1.)*rho)
-
-#--------------------------------------------------
-
-#File
-file = h5py.File(fileName, 'w')
-
-# Header
-grp = file.create_group("/Header")
-grp.attrs["BoxSize"] = boxSize
-grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
-grp.attrs["Time"] = 0.0
-grp.attrs["NumFilesPerSnapshot"] = 1
-grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-grp.attrs["Flag_Entropy_ICs"] = 0
-
-#Runtime parameters
-grp = file.create_group("/RuntimePars")
-grp.attrs["PeriodicBoundariesOn"] = periodic
-
-#Units
-grp = file.create_group("/Units")
-grp.attrs["Unit length in cgs (U_L)"] = 1.
-grp.attrs["Unit mass in cgs (U_M)"] = 1.
-grp.attrs["Unit time in cgs (U_t)"] = 1.
-grp.attrs["Unit current in cgs (U_I)"] = 1.
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.
-
-#Particle group
-grp = file.create_group("/PartType0")
-
-v = zeros((numPart, 3))
-ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
-ds[()] = v
-v = zeros(1)
-
-m = full((numPart, 1), mass)
-ds = grp.create_dataset('Masses', (numPart,1), 'f')
-ds[()] = m
-m = zeros(1)
-
-h = full((numPart, 1), eta * boxSize / L)
-ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
-ds[()] = h
-h = zeros(1)
-
-u = full((numPart, 1), internalEnergy)
-ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
-ds[()] = u
-u = zeros(1)
-
-
-ids = linspace(0, numPart, numPart, endpoint=False).reshape((numPart,1))
-ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
-ds[()] = ids + 1
-x = ids % L;
-y = ((ids - x) / L) % L;
-z = (ids - x - L * y) / L**2;
-coords = zeros((numPart, 3))
-coords[:,0] = z[:,0] * boxSize / L + boxSize / (2*L)
-coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L)
-coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L)
-ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
-ds[()] = coords
-
-file.close()
diff --git a/examples/IsothermalPotential/energy_plot.py b/examples/IsothermalPotential/energy_plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..0afa6fa93fa2a992e6ddeab3c9d33538c0b41de3
--- /dev/null
+++ b/examples/IsothermalPotential/energy_plot.py
@@ -0,0 +1,120 @@
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 12,
+'legend.fontsize': 12,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (3.15,3.15),
+'figure.subplot.left' : 0.145,
+'figure.subplot.right' : 0.99,
+'figure.subplot.bottom' : 0.11,
+'figure.subplot.top' : 0.99,
+'figure.subplot.wspace' : 0.15,
+'figure.subplot.hspace' : 0.12,
+'lines.markersize' : 6,
+'lines.linewidth' : 3.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+
+import numpy as np
+import h5py as h5
+import sys
+
+# File containing the total energy
+stats_filename = "./energy.txt"
+
+# First snapshot
+snap_filename = "Isothermal_000.hdf5"
+f = h5.File(snap_filename,'r')
+
+# Read the units parameters from the snapshot
+units = f["InternalCodeUnits"]
+unit_mass = units.attrs["Unit mass in cgs (U_M)"]
+unit_length = units.attrs["Unit length in cgs (U_L)"]
+unit_time = units.attrs["Unit time in cgs (U_t)"]
+
+# Read the header
+header = f["Header"]
+box_size = float(header.attrs["BoxSize"][0])
+
+# Read the properties of the potential
+parameters = f["Parameters"]
+R200 = 100
+Vrot = float(parameters.attrs["IsothermalPotential:vrot"])
+centre = [box_size/2, box_size/2, box_size/2]
+f.close()
+
+# Read the statistics summary
+file_energy = np.loadtxt("energy.txt")
+time_stats = file_energy[:,0]
+E_kin_stats = file_energy[:,3]
+E_pot_stats = file_energy[:,5]
+E_tot_stats = E_kin_stats + E_pot_stats
+
+# Read the snapshots
+time_snap = np.zeros(402)
+E_kin_snap = np.zeros(402)
+E_pot_snap = np.zeros(402)
+E_tot_snap = np.zeros(402)
+Lz_snap = np.zeros(402)
+
+# Read all the particles from the snapshots
+for i in range(402):
+ snap_filename = "Isothermal_%0.3d.hdf5"%i
+ f = h5.File(snap_filename,'r')
+
+ pos_x = f["PartType1/Coordinates"][:,0]
+ pos_y = f["PartType1/Coordinates"][:,1]
+ pos_z = f["PartType1/Coordinates"][:,2]
+ vel_x = f["PartType1/Velocities"][:,0]
+ vel_y = f["PartType1/Velocities"][:,1]
+ vel_z = f["PartType1/Velocities"][:,2]
+ mass = f["/PartType1/Masses"][:]
+
+ r = np.sqrt((pos_x[:] - centre[0])**2 + (pos_y[:] - centre[1])**2 + (pos_z[:] - centre[2])**2)
+ Lz = (pos_x[:] - centre[0]) * vel_y[:] - (pos_y[:] - centre[1]) * vel_x[:]
+
+ time_snap[i] = f["Header"].attrs["Time"]
+ E_kin_snap[i] = np.sum(0.5 * mass * (vel_x[:]**2 + vel_y[:]**2 + vel_z[:]**2))
+ E_pot_snap[i] = np.sum(-mass * Vrot**2 * log(r))
+ E_tot_snap[i] = E_kin_snap[i] + E_pot_snap[i]
+ Lz_snap[i] = np.sum(Lz)
+
+# Plot energy evolution
+figure()
+plot(time_stats, E_kin_stats, "r-", lw=0.5, label="Kinetic energy")
+plot(time_stats, E_pot_stats, "g-", lw=0.5, label="Potential energy")
+plot(time_stats, E_tot_stats, "k-", lw=0.5, label="Total energy")
+
+plot(time_snap[::10], E_kin_snap[::10], "rD", lw=0.5, ms=2)
+plot(time_snap[::10], E_pot_snap[::10], "gD", lw=0.5, ms=2)
+plot(time_snap[::10], E_tot_snap[::10], "kD", lw=0.5, ms=2)
+
+legend(loc="center right", fontsize=8, frameon=False, handlelength=3, ncol=1)
+xlabel("${\\rm{Time}}$", labelpad=0)
+ylabel("${\\rm{Energy}}$",labelpad=0)
+xlim(0, 8)
+
+savefig("energy.png", dpi=200)
+
+# Plot angular momentum evolution
+figure()
+plot(time_snap, Lz_snap, "k-", lw=0.5, ms=2)
+
+xlabel("${\\rm{Time}}$", labelpad=0)
+ylabel("${\\rm{Angular~momentum}}$",labelpad=0)
+xlim(0, 8)
+
+savefig("angular_momentum.png", dpi=200)
+
+
diff --git a/examples/IsothermalPotential/isothermal.yml b/examples/IsothermalPotential/isothermal.yml
index 0de99779f07591a5b71be11b75bc56ec741ddaed..8d9ec3875e405d95a89b3486bca5fd3465a3e20d 100644
--- a/examples/IsothermalPotential/isothermal.yml
+++ b/examples/IsothermalPotential/isothermal.yml
@@ -15,7 +15,7 @@ TimeIntegration:
# Parameters governing the conserved quantities statistics
Statistics:
- delta_time: 1e-2 # Time between statistics output
+ delta_time: 1e-3 # Time between statistics output
# Parameters governing the snapshots
Snapshots:
@@ -23,25 +23,18 @@ Snapshots:
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.02 # Time difference between consecutive outputs (in internal units)
-# Parameters for the hydrodynamics scheme
-SPH:
- resolution_eta: 1.2349 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
- delta_neighbours: 1. # The tolerance for the targetted number of neighbours.
- CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
- max_smoothing_length: 40. # Maximal smoothing length allowed (in internal units).
-
# Parameters related to the initial conditions
InitialConditions:
file_name: Isothermal.hdf5 # The file to read
- shift_x: 100. # A shift to apply to all particles read from the ICs (in internal units).
- shift_y: 100.
- shift_z: 100.
+ shift_x: 200. # Shift all particles to be in the potential
+ shift_y: 200.
+ shift_z: 200.
# External potential parameters
IsothermalPotential:
- position_x: 100. # location of centre of isothermal potential in internal units
- position_y: 100.
- position_z: 100.
+ position_x: 0. # location of centre of isothermal potential in internal units
+ position_y: 0.
+ position_z: 0.
vrot: 200. # rotation speed of isothermal potential in internal units
- timestep_mult: 0.03 # controls time step
-
+ timestep_mult: 0.01 # controls time step
+ epsilon: 0. # No softening at the centre of the halo
diff --git a/examples/IsothermalPotential/makeIC.py b/examples/IsothermalPotential/makeIC.py
index 976119f0a312c5acc81fab943ba3cf5769102269..7d1c5361f9a255365517226e49c55a8a50c4d6ce 100644
--- a/examples/IsothermalPotential/makeIC.py
+++ b/examples/IsothermalPotential/makeIC.py
@@ -30,10 +30,10 @@ import random
# all particles move in the xy plane, and start at y=0
# physical constants in cgs
-NEWTON_GRAVITY_CGS = 6.672e-8
+NEWTON_GRAVITY_CGS = 6.67408e-8
SOLAR_MASS_IN_CGS = 1.9885e33
PARSEC_IN_CGS = 3.0856776e18
-PROTON_MASS_IN_CGS = 1.6726231e24
+PROTON_MASS_IN_CGS = 1.672621898e24
YEAR_IN_CGS = 3.154e+7
# choice of units
@@ -66,17 +66,12 @@ N = int(sys.argv[1]) # Number of particles
icirc = int(sys.argv[2]) # if = 0, all particles are on circular orbits, if = 1, Lz/Lcirc uniform in ]0,1[
L = N**(1./3.)
-# these are not used but necessary for I/O
-rho = 2. # Density
-P = 1. # Pressure
-gamma = 5./3. # Gas adiabatic index
fileName = "Isothermal.hdf5"
#---------------------------------------------------
numPart = N
mass = 1
-internalEnergy = P / ((gamma - 1.)*rho)
#--------------------------------------------------
@@ -111,7 +106,6 @@ grp.attrs["PeriodicBoundariesOn"] = periodic
numpy.random.seed(1234)
#Particle group
-#grp0 = file.create_group("/PartType0")
grp1 = file.create_group("/PartType1")
#generate particle positions
radius = Radius * (numpy.random.rand(N))**(1./3.)
@@ -119,10 +113,8 @@ ctheta = -1. + 2 * numpy.random.rand(N)
stheta = numpy.sqrt(1.-ctheta**2)
phi = 2 * math.pi * numpy.random.rand(N)
r = numpy.zeros((numPart, 3))
-#r[:,0] = radius * stheta * numpy.cos(phi)
-#r[:,1] = radius * stheta * numpy.sin(phi)
-#r[:,2] = radius * ctheta
r[:,0] = radius
+
#
speed = vrot
v = numpy.zeros((numPart, 3))
@@ -146,17 +138,6 @@ ds = grp1.create_dataset('Masses', (numPart,), 'f')
ds[()] = m
m = numpy.zeros(1)
-h = numpy.full((numPart, ), 1.1255 * boxSize / L, dtype='f')
-ds = grp1.create_dataset('SmoothingLength', (numPart,), 'f')
-ds[()] = h
-h = numpy.zeros(1)
-
-u = numpy.full((numPart, ), internalEnergy, dtype='f')
-ds = grp1.create_dataset('InternalEnergy', (numPart,), 'f')
-ds[()] = u
-u = numpy.zeros(1)
-
-
ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False, dtype='L')
ds = grp1.create_dataset('ParticleIDs', (numPart, ), 'L')
ds[()] = ids
diff --git a/examples/IsothermalPotential/run.sh b/examples/IsothermalPotential/run.sh
index 28a3cc0910f986f84bcd603091543643356f1c4a..976fbddc01cf7a3dcbb114d437ddb8f03b4d54bd 100755
--- a/examples/IsothermalPotential/run.sh
+++ b/examples/IsothermalPotential/run.sh
@@ -7,4 +7,7 @@ then
python makeIC.py 1000 1
fi
+rm -rf Isothermal_*.hdf5
../swift -g -t 1 isothermal.yml 2>&1 | tee output.log
+
+python energy_plot.py
diff --git a/examples/IsothermalPotential/test.pro b/examples/IsothermalPotential/test.pro
deleted file mode 100644
index edfa50121d2e5adb7e039f3c38d6d4c0b4d5e34f..0000000000000000000000000000000000000000
--- a/examples/IsothermalPotential/test.pro
+++ /dev/null
@@ -1,168 +0,0 @@
-;
-; test energy / angular momentum conservation of test problem
-;
-
-iplot = 1 ; if iplot = 1, make plot of E/Lz conservation, else, simply compare final and initial energy
-
-; set physical constants
-@physunits
-
-indir = './'
-basefile = 'Isothermal_'
-
-; set properties of potential
-uL = 1e3 * phys.pc ; unit of length
-uM = phys.msun ; unit of mass
-uV = 1d5 ; unit of velocity
-vrot = 200. ; km/s
-r200 = 100. ; virial radius
-
-; derived units
-constG = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ;
-pcentre = [100.,100.,100.] * 1d3 * pc / uL
-
-;
-infile = indir + basefile + '*'
-spawn,'ls -1 '+infile,res
-nfiles = n_elements(res)
-
-
-
-; choose: calculate change of energy and Lz, comparing first and last
-; snapshots for all particles, or do so for a subset
-
-; compare all
-ifile = 0
-inf = indir + basefile + strtrim(string(ifile,'(i3.3)'),1) + '.hdf5'
-id = h5rd(inf,'PartType1/ParticleIDs')
-nfollow = n_elements(id)
-
-; follow a subset
-nfollow = 500 ; number of particles to follow
-
-;
-if (iplot eq 1) then begin
- nskip = 1
- nsave = nfiles
-endif else begin
- nskip = nfiles - 2
- nsave = 2
-endelse
-
-;
-lout = fltarr(nfollow, nsave) ; Lz
-xout = fltarr(nfollow, nsave) ; x
-yout = fltarr(nfollow, nsave) ; y
-zout = fltarr(nfollow, nsave) ; z
-eout = fltarr(nfollow, nsave) ; energies
-ekin = fltarr(nfollow, nsave)
-epot = fltarr(nfollow, nsave)
-tout = fltarr(nsave)
-
-
-
-ifile = 0
-isave = 0
-for ifile=0,nfiles-1,nskip do begin
- inf = indir + basefile + strtrim(string(ifile,'(i3.3)'),1) + '.hdf5'
- time = h5ra(inf, 'Header','Time')
- p = h5rd(inf,'PartType1/Coordinates')
- v = h5rd(inf,'PartType1/Velocities')
- id = h5rd(inf,'PartType1/ParticleIDs')
- indx = sort(id)
-;
- id = id[indx]
- for ic=0,2 do begin
- tmp = reform(p[ic,*]) & p[ic,*] = tmp[indx]
- tmp = reform(v[ic,*]) & v[ic,*] = tmp[indx]
- endfor
-
-
-; calculate energy
- dd = size(p,/dimen) & npart = dd[1]
- ener = fltarr(npart)
- dr = fltarr(npart) & dv = dr
- for ic=0,2 do dr[*] = dr[*] + (p[ic,*]-pcentre[ic])^2
- for ic=0,2 do dv[*] = dv[*] + v[ic,*]^2
- xout[*,isave] = p[0,0:nfollow-1]-pcentre[0]
- yout[*,isave] = p[1,0:nfollow-1]-pcentre[1]
- zout[*,isave] = p[2,0:nfollow-1]-pcentre[2]
- Lz = (p[0,*]-pcentre[0]) * v[1,*] - (p[1,*]-pcentre[1]) * v[0,*]
- dr = sqrt(dr)
-; print,'time = ',time,p[0,0],v[0,0],id[0]
- ek = 0.5 * dv
-; ep = - constG * mextern / dr
- ep = -vrot*vrot * (1 + alog(r200/dr))
- ener = ek + ep
- tout(isave) = time
- lout[*,isave] = lz[0:nfollow-1]
- eout(*,isave) = ener[0:nfollow-1]
- ekin(*,isave) = ek[0:nfollow-1]
- epot(*,isave) = ep[0:nfollow-1]
-
-; write some output
-; print,' time= ',time,' e= ',eout[0],' Lz= ',lz[0],format='(%a %f %a
-; %f)'
- print,format='('' time= '',f7.1,'' E= '',f9.2,'' Lz= '',e9.2)', time,eout[0],lz[0]
- isave = isave + 1
-
-endfor
-x0 = reform(xout[0,*])
-y0 = reform(xout[1,*])
-z0 = reform(xout[2,*])
-
-; calculate relative energy change
-de = 0.0 * eout
-dl = 0.0 * lout
-nsave = isave
-for ifile=1, nsave-1 do de[*,ifile] = (eout[*,ifile]-eout[*,0])/eout[*,0]
-for ifile=1, nsave-1 do dl[*,ifile] = (lout[*,ifile] - lout[*,0])/lout[*,0]
-
-
-; calculate statistics of energy changes
-print,' relatve energy change: (per cent) ',minmax(de) * 100.
-print,' relative Lz change: (per cent) ',minmax(dl) * 100.
-
-; plot enery and Lz conservation for some particles
-if(iplot eq 1) then begin
-; plot results on energy conservation for some particles
- nplot = min(10, nfollow)
- win,0
- xr = [min(tout), max(tout)]
- yr = [-2,2]*1d-2 ; in percent
- plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/nodata,xtitle='time',ytitle='dE/E, dL/L (%)'
- for i=0,nplot-1 do oplot,tout,de[i,*]
- for i=0,nplot-1 do oplot,tout,dl[i,*],color=red
- legend,['dE/E','dL/L'],linestyle=[0,0],color=[black,red],box=0,/bottom,/left
- screen_to_png,'e-time.png'
-
-; plot orbits of those particles
- win,2
- xr = [-100,100]
- yr = xr
- plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/iso,/nodata,xtitle='x',ytitle='y'
- color = floor(findgen(nplot)*255/float(nplot))
- for i=0,nplot-1 do oplot,xout[i,*],yout[i,*],color=color(i)
- screen_to_png,'orbit.png'
-
-; plot radial position of these particles
- win,4
- xr = [min(tout), max(tout)]
- yr = [0,80]
- plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/nodata,xtitle='t',ytitle='r'
- color = floor(findgen(nplot)*255/float(nplot))
-for i=0,nplot-1 do begin dr = sqrt(reform(xout[i,*])^2 + reform(yout[i,*])^2) & oplot,tout,dr,color=color[i] & endfor
- screen_to_png,'r-time.png'
-
-; make histogram of energy changes at end
- win,6
- ohist,de,x,y,-0.05,0.05,0.001
- plot,x,y,psym=10,xtitle='de (%)'
- screen_to_png,'de-hist.png'
-
-
-endif
-
-end
-
-