Merge branch 'master' into Self_gravity_BH

.gitignore

@@ -28,10 +28,13 @@ examples/*.xmf
 examples/used_parameters.yml
 examples/energy.txt
 examples/*/*.xmf
-examples/*/*.dat
 examples/*/*.hdf5
+examples/*/*.txt
 examples/*/used_parameters.yml
-examples/*/energy.txt
+examples/*/*/*.xmf
+examples/*/*/*.hdf5
+examples/*/*/*.txt
+examples/*/*/used_parameters.yml
 
 tests/testPair
 tests/brute_force_standard.dat
@@ -39,10 +42,13 @@ tests/swift_dopair_standard.dat
 tests/brute_force_perturbed.dat
 tests/swift_dopair_perturbed.dat
 tests/test27cells
+tests/test125cells
 tests/brute_force_27_standard.dat
 tests/swift_dopair_27_standard.dat
 tests/brute_force_27_perturbed.dat
 tests/swift_dopair_27_perturbed.dat
+tests/brute_force_125_standard.dat
+tests/swift_dopair_125_standard.dat
 tests/testGreetings
 tests/testReading
 tests/input.hdf5
@@ -52,10 +58,14 @@ tests/testSPHStep
 tests/testKernel
 tests/testKernelGrav
 tests/testFFT
+tests/testInteractions
+tests/testSymmetry
+tests/testMaths
 tests/testParser
 tests/parser_output.yml
 tests/test27cells.sh
 tests/test27cellsPerturbed.sh
+tests/test125cells.sh
 tests/testPair.sh
 tests/testPairPerturbed.sh
 tests/testParser.sh

configure.ac

@@ -487,6 +487,7 @@ AC_CONFIG_FILES([tests/testPair.sh], [chmod +x tests/testPair.sh])
 AC_CONFIG_FILES([tests/testPairPerturbed.sh], [chmod +x tests/testPairPerturbed.sh])
 AC_CONFIG_FILES([tests/test27cells.sh], [chmod +x tests/test27cells.sh])
 AC_CONFIG_FILES([tests/test27cellsPerturbed.sh], [chmod +x tests/test27cellsPerturbed.sh])
+AC_CONFIG_FILES([tests/test125cells.sh], [chmod +x tests/test125cells.sh])
 AC_CONFIG_FILES([tests/testParser.sh], [chmod +x tests/testParser.sh])
 
 # Report general configuration.

examples/CosmoVolume/cosmoVolume.yml

@@ -23,11 +23,6 @@ Snapshots:
   basename:            cosmo # Common part of the name of output files
   time_first:          0.    # Time of the first output (in internal units)
   delta_time:          0.05  # Time difference between consecutive outputs (in internal units)
-  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 conserved quantities statistics
 Statistics:

examples/EAGLE_12/README

@@ -11,3 +11,6 @@ integration.
 
 The particle load of the main EAGLE simulation can be reproduced by
 running these ICs on 8 cores.
+
+MD5 checksum of the ICs: 
+88877f5bb0ee21488c20b8f065fc74db  EAGLE_ICs_12.hdf5

examples/EAGLE_12/eagle_12.yml

 # 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
+  UnitMass_in_cgs:     1.989e43      # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.085678e24   # Mpc in centimeters
+  UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
+  UnitCurrent_in_cgs:  1             # Amperes
+  UnitTemp_in_cgs:     1             # Kelvin
 
 # Parameters for the task scheduling
 Scheduler:
@@ -14,20 +14,15 @@ Scheduler:
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
-  time_end:   1.    # 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-2  # The maximal time-step size of the simulation (in internal units).
+  time_end:   1e-2  # The end time of the simulation (in internal units).
+  dt_min:     1e-10 # 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:            eagle # Common part of the name of output files
   time_first:          0.    # Time of the first output (in internal units)
-  delta_time:          0.05  # Time difference between consecutive outputs (in internal units)
-  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
+  delta_time:          1e-3  # Time difference between consecutive outputs (in internal units)
 
 # Parameters governing the conserved quantities statistics
 Statistics:

examples/EAGLE_25/README

@@ -11,3 +11,6 @@ Everything is ready to be run without cosmological integration.
 
 The particle load of the main EAGLE simulation can be reproduced by
 running these ICs on 64 cores.
+
+MD5 checksum of the ICs: 
+ada2c728db2bd2d77a20c4eef52dfaf1  EAGLE_ICs_25.hdf5

examples/EAGLE_25/eagle_25.yml

 # 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
+  UnitMass_in_cgs:     1.989e43      # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.085678e24   # Mpc in centimeters
+  UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
+  UnitCurrent_in_cgs:  1             # Amperes
+  UnitTemp_in_cgs:     1             # Kelvin
 
 # Parameters for the task scheduling
 Scheduler:
@@ -14,20 +14,15 @@ Scheduler:
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
-  time_end:   1.    # 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-2  # The maximal time-step size of the simulation (in internal units).
+  time_end:   1e-2  # The end time of the simulation (in internal units).
+  dt_min:     1e-10 # 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:            eagle # Common part of the name of output files
   time_first:          0.    # Time of the first output (in internal units)
-  delta_time:          0.05  # Time difference between consecutive outputs (in internal units)
-  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
+  delta_time:          1e-3  # Time difference between consecutive outputs (in internal units)
 
 # Parameters governing the conserved quantities statistics
 Statistics:

examples/EAGLE_50/README

@@ -11,3 +11,6 @@ Everything is ready to be run without cosmological integration.
 
 The particle load of the main EAGLE simulation can be reproduced by
 running these ICs on 512 cores.
+
+MD5 checksum of the ICs: 
+3591b579bd108ddf0e555092bdfbf97f  EAGLE_ICs_50.hdf5

examples/EAGLE_50/eagle_50.yml

 # 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
+  UnitMass_in_cgs:     1.989e43      # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.085678e24   # Mpc in centimeters
+  UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
+  UnitCurrent_in_cgs:  1             # Amperes
+  UnitTemp_in_cgs:     1             # Kelvin
 
 # Parameters for the task scheduling
 Scheduler:
@@ -14,20 +14,15 @@ Scheduler:
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
-  time_end:   1.    # 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-2  # The maximal time-step size of the simulation (in internal units).
+  time_end:   1e-2  # The end time of the simulation (in internal units).
+  dt_min:     1e-10 # 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:            eagle # Common part of the name of output files
   time_first:          0.    # Time of the first output (in internal units)
-  delta_time:          0.05  # Time difference between consecutive outputs (in internal units)
-  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
+  delta_time:          1e-3  # Time difference between consecutive outputs (in internal units)
 
 # Parameters governing the conserved quantities statistics
 Statistics:

examples/ExternalPointMass/externalPointMass.yml

@@ -18,11 +18,6 @@ Snapshots:
   basename:            pointMass # Common part of the name of output files
   time_first:          0.        # Time of the first output (in internal units)
   delta_time:          0.02      # Time difference between consecutive outputs (in internal units)
-  UnitMass_in_cgs:     1.9885e33     # Grams
-  UnitLength_in_cgs:   3.0856776e21  # Centimeters
-  UnitVelocity_in_cgs: 1e5           # Centimeters per second
-  UnitCurrent_in_cgs:  1   # Amperes
-  UnitTemp_in_cgs:     1   # Kelvin
 
 # Parameters governing the conserved quantities statistics
 Statistics:
@@ -48,4 +43,5 @@ PointMass:
   position_y:      50.
   position_z:      50.	
   mass:            1e10     # mass of external point mass in internal units
+  timestep_mult:   0.03     # controls time step
 

examples/IsothermalPotential/GravityOnly/README

+;
+; this probelm generates a set of gravity particles in an isothermal
+; potential and follows their orbits. Tests verify consdevation of
+; energy and angular momentum
+;
+;

examples/IsothermalPotential/GravityOnly/isothermal.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.9885e33     # Grams
+  UnitLength_in_cgs:   3.0856776e21  # Centimeters
+  UnitVelocity_in_cgs: 1e5           # 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:   8.    # 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-1  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-2 # Time between statistics output
+  
+# Parameters governing the snapshots
+Snapshots:
+  basename:            Isothermal # Common part of the name of output files
+  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.
+ 
+# External potential parameters
+IsothermalPotential:
+  position_x:      100.     # location of centre of isothermal potential in internal units
+  position_y:      100.
+  position_z:      100.	
+  vrot:            200.     # rotation speed of isothermal potential in internal units
+  timestep_mult:   0.03     # controls time step
+

examples/IsothermalPotential/GravityOnly/makeIC.py

+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ #                    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
+import numpy
+import math
+import random
+
+# Generates N particles in a spherical distribution centred on [0,0,0], to be moved in an isothermal potential
+# usage: python makeIC.py 1000 0 : generate 1000 particles on circular orbits
+#        python makeIC.py 1000 1 : generate 1000 particles with Lz/L uniform in [0,1]
+# all particles move in the xy plane, and start at y=0
+
+# physical constants in cgs
+NEWTON_GRAVITY_CGS  = 6.672e-8
+SOLAR_MASS_IN_CGS   = 1.9885e33
+PARSEC_IN_CGS       = 3.0856776e18
+PROTON_MASS_IN_CGS  = 1.6726231e24
+YEAR_IN_CGS         = 3.154e+7
+
+# choice of units
+const_unit_length_in_cgs   =   (1000*PARSEC_IN_CGS)
+const_unit_mass_in_cgs     =   (SOLAR_MASS_IN_CGS)
+const_unit_velocity_in_cgs =   (1e5)
+
+print "UnitMass_in_cgs:     ", const_unit_mass_in_cgs 
+print "UnitLength_in_cgs:   ", const_unit_length_in_cgs
+print "UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs
+
+
+# rotation speed of isothermal potential [km/s]
+vrot_kms = 200.
+
+# derived units
+const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs)
+const_G                = ((NEWTON_GRAVITY_CGS*const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/(const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs)))
+print 'G=', const_G
+vrot   = vrot_kms * 1e5 / const_unit_velocity_in_cgs
+
+
+# Parameters
+periodic= 1            # 1 For periodic box
+boxSize = 400.         #  [kpc]
+Radius  = 100.          # maximum radius of particles [kpc]
+G       = const_G 
+
+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)
+
+#--------------------------------------------------
+
+#File
+file = h5py.File(fileName, 'w')
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = const_unit_length_in_cgs
+grp.attrs["Unit mass in cgs (U_M)"] = const_unit_mass_in_cgs 
+grp.attrs["Unit time in cgs (U_t)"] = const_unit_length_in_cgs / const_unit_velocity_in_cgs
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+grp.attrs["NumPart_Total"] =  [0, numPart, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [0, numPart, 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, 0, 0, 0, 0, 0]
+
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+# set seed for random number
+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.) 
+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))
+omega  = speed / radius
+period = 2.*math.pi/omega
+print 'period = minimum = ',min(period), ' maximum = ',max(period)
+
+omegav = omega
+if (icirc != 0):
+    omegav = omega * numpy.random.rand(N)
+
+v[:,0] = -omegav * r[:,1]
+v[:,1] =  omegav * r[:,0]
+
+ds = grp1.create_dataset('Velocities', (numPart, 3), 'f')
+ds[()] = v
+v = numpy.zeros(1)
+
+m = numpy.full((numPart, ), mass)
+ds = grp1.create_dataset('Masses', (numPart,), 'f')
+ds[()] = m
+m = numpy.zeros(1)
+
+h = numpy.full((numPart, ), 1.1255 * boxSize / L)
+ds = grp1.create_dataset('SmoothingLength', (numPart,), 'f')
+ds[()] = h
+h = numpy.zeros(1)
+
+u = numpy.full((numPart, ), internalEnergy)
+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
+
+ds = grp1.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = r
+
+
+file.close()

examples/IsothermalPotential/GravityOnly/run.sh

+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e Isothermal.hdf5 ]
+then
+    echo "Generating initial conditions for the isothermal potential box example..."
+    python makeIC.py 1000 1
+fi
+
+../../swift -g -t 2 isothermal.yml

examples/IsothermalPotential/GravityOnly/test.pro

+;
+;  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