diff --git a/examples/Disk-Patch/HydroStatic/disk-patch-icc.yml b/examples/Disk-Patch/HydroStatic/disk-patch-icc.yml
new file mode 100644
index 0000000000000000000000000000000000000000..584a8df9c297dbb1337a2cb343f5eb5d9fcdbeba
--- /dev/null
+++ b/examples/Disk-Patch/HydroStatic/disk-patch-icc.yml
@@ -0,0 +1,61 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.9885e33 # Grams
+ UnitLength_in_cgs: 3.0856776e18 # 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: 968. # The end time of the simulation (in internal units).
+ dt_min: 1e-4 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1. # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1 # Time between statistics output
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: Disk-Patch # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 12. # 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: 0.1 # The tolerance for the targetted number of neighbours.
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ max_ghost_iterations: 30 # Maximal number of iterations allowed to converge towards the smoothing length.
+ max_smoothing_length: 70. # Maximal smoothing length allowed (in internal units).
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: Disk-Patch.hdf5 # The file to read
+ h_scaling: 1. # A scaling factor to apply to all smoothing lengths in the ICs.
+ shift_x: 0. # A shift to apply to all particles read from the ICs (in internal units).
+ shift_y: 0.
+ shift_z: 0.
+
+# External potential parameters
+PointMass:
+ position_x: 50. # location of external point mass in internal units
+ position_y: 50.
+ position_z: 50.
+ mass: 1e10 # mass of external point mass in internal units
+ timestep_mult: 0.03 # controls time step
+
+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
+
+Disk-PatchPotential:
+ surface_density: 10.
+ scale_height: 100.
+ z_disk: 200.
+ timestep_mult: 0.03
diff --git a/examples/Disk-Patch/HydroStatic/dynamic.pro b/examples/Disk-Patch/HydroStatic/dynamic.pro
index b598530d766d5c092802b9b0acc0aae9d067d8c6..c02c65fe418e84cdd62978dbddcf5a641fa4c156 100644
--- a/examples/Disk-Patch/HydroStatic/dynamic.pro
+++ b/examples/Disk-Patch/HydroStatic/dynamic.pro
@@ -19,9 +19,14 @@ uV = 1d5 ; unit of velocity
surface_density = 10.
scale_height = 100.
z_disk = 200.;
+gamma = 5./3.
+
; derived units
constG = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ;
pcentre = [0.,0.,z_disk] * pc / uL
+utherm = !pi * constG * surface_density * scale_height / (gamma-1.)
+soundspeed = sqrt(gamma * (gamma-1.) * utherm)
+t_dyn = sqrt(scale_height / (constG * surface_density))
;
infile = indir + basefile + '*'
@@ -48,7 +53,7 @@ rbins = (surface_density/(2.*scale_height)) / cosh(abs(zbins)/scale_height)^2
; plot analytic profile
wset,0
plot,[0],[0],xr=[0,2*scale_height],yr=[0,max(rbins)],/nodata,xtitle='|z|',ytitle=textoidl('\rho')
-oplot,zbins,rbins
+oplot,zbins,rbins,color=blue
ifile = 0
nskip = nfiles - 1
@@ -56,8 +61,10 @@ isave = 0
nplot = 8192 ; randomly plot particles
color = floor(findgen(nfiles)/float(nfiles-1)*255)
;for ifile=0,nfiles-1,nskip do begin
-tsave = [0.]
-for ifile=0,nfiles-1,nskip do begin
+tsave = [0.]
+toplot = [1,nfiles-1]
+for iplot=0,1 do begin
+ ifile = toplot[iplot]
inf = indir + basefile + strtrim(string(ifile,'(i3.3)'),1) + '.hdf5'
time = h5ra(inf, 'Header','Time')
tsave = [tsave, time]
@@ -86,19 +93,37 @@ for ifile=0,nfiles-1,nskip do begin
ip = floor(randomu(ifile+1,nplot)*n_elements(rho))
color = red
- if(ifile eq 0) then color=black
+ if(ifile eq 1) then begin
+ color=black
+ endif else begin
+ color=red
+ endelse
oplot,abs(p[2,ip]), rho[ip], psym=3, color=color
isave = isave + 1
endfor
-tsave = tsave[1:*]
+
+; time in units of dynamical time
+tsave = tsave[1:*] / t_dyn
label = ['']
-for i=0,n_elements(tsave)-1 do label=[label,'t='+string(tsave[i],format='(f8.0)')]
+for i=0,n_elements(tsave)-1 do label=[label,'time/t_dynamic='+string(tsave[i],format='(f8.0)')]
label = label[1:*]
-legend,[label[0],label[1]],linestyle=[0,0],color=[black,red],box=0,/top,/right
-
+legend,['analytic',label[0],label[1]],linestyle=[0,0,0],color=[blue,black,red],box=0,/top,/right
+
+; make histograms of particle velocities
+xr = 1d-3 * [-1,1]
+bsize = 1.d-5
+ohist,v[0,*]/soundspeed,x,vx,xr[0],xr[1],bsize
+ohist,v[1,*]/soundspeed,y,vy,xr[0],xr[1],bsize
+ohist,v[2,*]/soundspeed,z,vz,xr[0],xr[1],bsize
+wset,2
+plot,x,vx,psym=10,xtitle='velocity/soundspeed',ytitle='pdf',/nodata,xr=xr,/xs
+oplot,x,vx,psym=10,color=black
+oplot,y,vy,psym=10,color=blue
+oplot,z,vz,psym=10,color=red
+legend,['vx/c','vy/c','vz/c'],linestyle=[0,0,0],color=[black,blue,red],box=0,/top,/right
end