diff --git a/examples/Disk-Patch/GravityOnly/disk-patch.yml b/examples/Disk-Patch/GravityOnly/disk-patch.yml
new file mode 100644
index 0000000000000000000000000000000000000000..37752a7a82e82c5d86b2603e3a30128c4d2664c6
--- /dev/null
+++ b/examples/Disk-Patch/GravityOnly/disk-patch.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: 480. # The end time of the simulation (in internal units).
+ dt_min: 1e-3 # 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.0 # 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: 8. # 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_ghost_iterations: 30 # Maximal number of iterations allowed to converge towards the smoothing length.
+ max_smoothing_length: 40. # 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: 300.
+ timestep_mult: 0.03
diff --git a/examples/Disk-Patch/GravityOnly/legend.pro b/examples/Disk-Patch/GravityOnly/legend.pro
new file mode 100755
index 0000000000000000000000000000000000000000..e510bdfcd7dcabed796cb090c6fa7d54536608b6
--- /dev/null
+++ b/examples/Disk-Patch/GravityOnly/legend.pro
@@ -0,0 +1,459 @@
+;+
+; NAME:
+; LEGEND
+; PURPOSE:
+; Create an annotation legend for a plot.
+; EXPLANATION:
+; This procedure makes a legend for a plot. The legend can contain
+; a mixture of symbols, linestyles, Hershey characters (vectorfont),
+; and filled polygons (usersym). A test procedure, legendtest.pro,
+; shows legend's capabilities. Placement of the legend is controlled
+; with keywords like /right, /top, and /center or by using a position
+; keyword for exact placement (position=[x,y]) or via mouse (/position).
+; CALLING SEQUENCE:
+; LEGEND [,items][,keyword options]
+; EXAMPLES:
+; The call:
+; legend,['Plus sign','Asterisk','Period'],psym=[1,2,3]
+; produces:
+; -----------------
+; | |
+; | + Plus sign |
+; | * Asterisk |
+; | . Period |
+; | |
+; -----------------
+; Each symbol is drawn with a plots command, so they look OK.
+; Other examples are given in optional output keywords.
+;
+; lines = indgen(6) ; for line styles
+; items = 'linestyle '+strtrim(lines,2) ; annotations
+; legend,items,linestyle=lines ; vertical legend---upper left
+; items = ['Plus sign','Asterisk','Period']
+; sym = [1,2,3]
+; legend,items,psym=sym ; ditto except using symbols
+; legend,items,psym=sym,/horizontal ; horizontal format
+; legend,items,psym=sym,box=0 ; sans border
+; legend,items,psym=sym,delimiter='=' ; embed '=' betw psym & text
+; legend,items,psym=sym,margin=2 ; 2-character margin
+; legend,items,psym=sym,position=[x,y] ; upper left in data coords
+; legend,items,psym=sym,pos=[x,y],/norm ; upper left in normal coords
+; legend,items,psym=sym,pos=[x,y],/device ; upper left in device coords
+; legend,items,psym=sym,/position ; interactive position
+; legend,items,psym=sym,/right ; at upper right
+; legend,items,psym=sym,/bottom ; at lower left
+; legend,items,psym=sym,/center ; approximately near center
+; legend,items,psym=sym,number=2 ; plot two symbols, not one
+; legend,items,/fill,psym=[8,8,8],colors=[10,20,30]; 3 filled squares
+; INPUTS:
+; items = text for the items in the legend, a string array.
+; For example, items = ['diamond','asterisk','square'].
+; You can omit items if you don't want any text labels.
+; OPTIONAL INPUT KEYWORDS:
+;
+; linestyle = array of linestyle numbers If linestyle[i] < 0, then omit
+; ith symbol or line to allow a multi-line entry. If
+; linestyle = -99 then text will be left-justified.
+; psym = array of plot symbol numbers. If psym[i] is negative, then a
+; line connects pts for ith item. If psym[i] = 8, then the
+; procedure usersym is called with vertices define in the
+; keyword usersym. If psym[i] = 88, then use the previously
+; defined user symbol
+; vectorfont = vector-drawn characters for the sym/line column, e.g.,
+; ['!9B!3','!9C!3','!9D!3'] produces an open square, a checkmark,
+; and a partial derivative, which might have accompanying items
+; ['BOX','CHECK','PARTIAL DERIVATIVE'].
+; There is no check that !p.font is set properly, e.g., -1 for
+; X and 0 for PostScript. This can produce an error, e.g., use
+; !20 with PostScript and !p.font=0, but allows use of Hershey
+; *AND* PostScript fonts together.
+; N. B.: Choose any of linestyle, psym, and/or vectorfont. If none is
+; present, only the text is output. If more than one
+; is present, all need the same number of elements, and normal
+; plot behaviour occurs.
+; By default, if psym is positive, you get one point so there is
+; no connecting line. If vectorfont[i] = '',
+; then plots is called to make a symbol or a line, but if
+; vectorfont[i] is a non-null string, then xyouts is called.
+; /help = flag to print header
+; /horizontal = flag to make the legend horizontal
+; /vertical = flag to make the legend vertical (D=vertical)
+; box = flag to include/omit box around the legend (D=include)
+; clear = flag to clear the box area before drawing the legend
+; delimiter = embedded character(s) between symbol and text (D=none)
+; colors = array of colors for plot symbols/lines (D=!P.color)
+; textcolors = array of colors for text (D=!P.color)
+; margin = margin around text measured in characters and lines
+; spacing = line spacing (D=bit more than character height)
+; pspacing = psym spacing (D=3 characters) (when number of symbols is
+; greater than 1)
+; charsize = just like !p.charsize for plot labels
+; charthick = just like !p.charthick for plot labels
+; thick = array of line thickness numbers (D = !P.thick), if used, then
+; linestyle must also be specified
+; position = data coordinates of the /top (D) /left (D) of the legend
+; normal = use normal coordinates for position, not data
+; device = use device coordinates for position, not data
+; number = number of plot symbols to plot or length of line (D=1)
+; usersym = 2-D array of vertices, cf. usersym in IDL manual.
+; (/USERSYM =square, default is to use existing USERSYM definition)
+; /fill = flag to fill the usersym
+; /left_legend = flag to place legend snug against left side of plot
+; window (D)
+; /right_legend = flag to place legend snug against right side of plot
+; window. If /right,pos=[x,y], then x is position of RHS and
+; text runs right-to-left.
+; /top_legend = flag to place legend snug against top of plot window (D)
+; /bottom = flag to place legend snug against bottom of plot window
+; /top,pos=[x,y] and /bottom,pos=[x,y] produce same positions.
+;
+; If LINESTYLE, PSYM, VECTORFONT, THICK, COLORS, or TEXTCOLORS are
+; supplied as scalars, then the scalar value is set for every line or
+; symbol in the legend.
+; Outputs:
+; legend to current plot device
+; OPTIONAL OUTPUT KEYWORDS:
+; corners = 4-element array, like !p.position, of the normalized
+; coords for the box (even if box=0): [llx,lly,urx,ury].
+; Useful for multi-column or multi-line legends, for example,
+; to make a 2-column legend, you might do the following:
+; c1_items = ['diamond','asterisk','square']
+; c1_psym = [4,2,6]
+; c2_items = ['solid','dashed','dotted']
+; c2_line = [0,2,1]
+; legend,c1_items,psym=c1_psym,corners=c1,box=0
+; legend,c2_items,line=c2_line,corners=c2,box=0,pos=[c1[2],c1[3]]
+; c = [c1[0]<c2[0],c1[1]<c2[1],c1[2]>c2[2],c1[3]>c2[3]]
+; plots,[c[0],c[0],c[2],c[2],c[0]],[c[1],c[3],c[3],c[1],c[1]],/norm
+; Useful also to place the legend. Here's an automatic way to place
+; the legend in the lower right corner. The difficulty is that the
+; legend's width is unknown until it is plotted. In this example,
+; the legend is plotted twice: the first time in the upper left, the
+; second time in the lower right.
+; legend,['1','22','333','4444'],linestyle=indgen(4),corners=corners
+; ; BOGUS LEGEND---FIRST TIME TO REPORT CORNERS
+; xydims = [corners[2]-corners[0],corners[3]-corners[1]]
+; ; SAVE WIDTH AND HEIGHT
+; chdim=[!d.x_ch_size/float(!d.x_size),!d.y_ch_size/float(!d.y_size)]
+; ; DIMENSIONS OF ONE CHARACTER IN NORMALIZED COORDS
+; pos = [!x.window[1]-chdim[0]-xydims[0] $
+; ,!y.window[0]+chdim[1]+xydims[1]]
+; ; CALCULATE POSITION FOR LOWER RIGHT
+; plot,findgen(10) ; SIMPLE PLOT; YOU DO WHATEVER YOU WANT HERE.
+; legend,['1','22','333','4444'],linestyle=indgen(4),pos=pos
+; ; REDO THE LEGEND IN LOWER RIGHT CORNER
+; You can modify the pos calculation to place the legend where you
+; want. For example to place it in the upper right:
+; pos = [!x.window[1]-chdim[0]-xydims[0],!y.window[1]-xydims[1]]
+; Common blocks:
+; none
+; Procedure:
+; If keyword help is set, call doc_library to print header.
+; See notes in the code. Much of the code deals with placement of the
+; legend. The main problem with placement is not being
+; able to sense the length of a string before it is output. Some crude
+; approximations are used for centering.
+; Restrictions:
+; Here are some things that aren't implemented.
+; - An orientation keyword would allow lines at angles in the legend.
+; - An array of usersyms would be nice---simple change.
+; - An order option to interchange symbols and text might be nice.
+; - Somebody might like double boxes, e.g., with box = 2.
+; - Another feature might be a continuous bar with ticks and text.
+; - There are no guards to avoid writing outside the plot area.
+; - There is no provision for multi-line text, e.g., '1st line!c2nd line'
+; Sensing !c would be easy, but !c isn't implemented for PostScript.
+; A better way might be to simply output the 2nd line as another item
+; but without any accompanying symbol or linestyle. A flag to omit
+; the symbol and linestyle is linestyle[i] = -1.
+; - There is no ability to make a title line containing any of titles
+; for the legend, for the symbols, or for the text.
+; Side Effects:
+; Modification history:
+; write, 24-25 Aug 92, F K Knight (knight@ll.mit.edu)
+; allow omission of items or omission of both psym and linestyle, add
+; corners keyword to facilitate multi-column legends, improve place-
+; ment of symbols and text, add guards for unequal size, 26 Aug 92, FKK
+; add linestyle(i)=-1 to suppress a single symbol/line, 27 Aug 92, FKK
+; add keyword vectorfont to allow characters in the sym/line column,
+; 28 Aug 92, FKK
+; add /top, /bottom, /left, /right keywords for automatic placement at
+; the four corners of the plot window. The /right keyword forces
+; right-to-left printing of menu. 18 Jun 93, FKK
+; change default position to data coords and add normal, data, and
+; device keywords, 17 Jan 94, FKK
+; add /center keyword for positioning, but it is not precise because
+; text string lengths cannot be known in advance, 17 Jan 94, FKK
+; add interactive positioning with /position keyword, 17 Jan 94, FKK
+; allow a legend with just text, no plotting symbols. This helps in
+; simply describing a plot or writing assumptions done, 4 Feb 94, FKK
+; added thick, symsize, and clear keyword Feb 96, W. Landsman HSTX
+; David Seed, HR Wallingford, d.seed@hrwallingford.co.uk
+; allow scalar specification of keywords, Mar 96, W. Landsman HSTX
+; added charthick keyword, June 96, W. Landsman HSTX
+; Made keyword names left,right,top,bottom,center longer,
+; Aug 16, 2000, Kim Tolbert
+; Added ability to have regular text lines in addition to plot legend
+; lines in legend. If linestyle is -99 that item is left-justified.
+; Previously, only option for no sym/line was linestyle=-1, but then text
+; was lined up after sym/line column. 10 Oct 2000, Kim Tolbert
+; Make default value of thick = !P.thick W. Landsman Jan. 2001
+; Don't overwrite existing USERSYM definition W. Landsman Mar. 2002
+;-
+pro legend, items, BOTTOM_LEGEND=bottom, BOX = box, CENTER_LEGEND=center, $
+ CHARTHICK=charthick, CHARSIZE = charsize, CLEAR = clear, COLORS = colorsi, $
+ CORNERS = corners, DATA=data, DELIMITER=delimiter, DEVICE=device, $
+ FILL=fill, HELP = help, HORIZONTAL=horizontal,LEFT_LEGEND=left, $
+ LINESTYLE=linestylei, MARGIN=margin, NORMAL=normal, NUMBER=number, $
+ POSITION=position,PSPACING=pspacing, PSYM=psymi, RIGHT_LEGEND=right, $
+ SPACING=spacing, SYMSIZE=symsize, TEXTCOLORS=textcolorsi, THICK=thicki, $
+ TOP_LEGEND=top, USERSYM=usersym, VECTORFONT=vectorfonti, VERTICAL=vertical
+;
+; =====>> HELP
+;
+on_error,2
+if keyword_set(help) then begin & doc_library,'legend' & return & endif
+;
+; =====>> SET DEFAULTS FOR SYMBOLS, LINESTYLES, AND ITEMS.
+;
+ ni = n_elements(items)
+ np = n_elements(psymi)
+ nl = n_elements(linestylei)
+ nth = n_elements(thicki)
+ nv = n_elements(vectorfonti)
+ nlpv = max([np,nl,nv])
+ n = max([ni,np,nl,nv]) ; NUMBER OF ENTRIES
+strn = strtrim(n,2) ; FOR ERROR MESSAGES
+if n eq 0 then message,'No inputs! For help, type legend,/help.'
+if ni eq 0 then begin
+ items = replicate('',n) ; DEFAULT BLANK ARRAY
+endif else begin
+ if size(items,/TNAME) NE 'STRING' then message, $
+ 'First parameter must be a string array. For help, type legend,/help.'
+ if ni ne n then message,'Must have number of items equal to '+strn
+endelse
+symline = (np ne 0) or (nl ne 0) ; FLAG TO PLOT SYM/LINE
+ if (np ne 0) and (np ne n) and (np NE 1) then message, $
+ 'Must have 0, 1 or '+strn+' elements in PSYM array.'
+ if (nl ne 0) and (nl ne n) and (nl NE 1) then message, $
+ 'Must have 0, 1 or '+strn+' elements in LINESTYLE array.'
+ if (nth ne 0) and (nth ne n) and (nth NE 1) then message, $
+ 'Must have 0, 1 or '+strn+' elements in THICK array.'
+
+ case nl of
+ 0: linestyle = intarr(n) ;Default = solid
+ 1: linestyle = intarr(n) + linestylei
+ else: linestyle = linestylei
+ endcase
+
+ case nth of
+ 0: thick = replicate(!p.thick,n) ;Default = !P.THICK
+ 1: thick = intarr(n) + thicki
+ else: thick = thicki
+ endcase
+
+ case np of ;Get symbols
+ 0: psym = intarr(n) ;Default = solid
+ 1: psym = intarr(n) + psymi
+ else: psym = psymi
+ endcase
+
+ case nv of
+ 0: vectorfont = replicate('',n)
+ 1: vectorfont = replicate(vectorfonti,n)
+ else: vectorfont = vectorfonti
+ endcase
+;
+; =====>> CHOOSE VERTICAL OR HORIZONTAL ORIENTATION.
+;
+if n_elements(horizontal) eq 0 then begin ; D=VERTICAL
+ if n_elements(vertical) eq 0 then vertical = 1
+endif else begin
+ if n_elements(vertical) eq 0 then vertical = not horizontal
+endelse
+;
+; =====>> SET DEFAULTS FOR OTHER OPTIONS.
+;
+if n_elements(box) eq 0 then box = 1
+if n_elements(clear) eq 0 then clear = 0
+
+if n_elements(margin) eq 0 then margin = 0.5
+if n_elements(delimiter) eq 0 then delimiter = ''
+if n_elements(charsize) eq 0 then charsize = !p.charsize
+if n_elements(charthick) eq 0 then charthick = !p.charthick
+if charsize eq 0 then charsize = 1
+if (n_elements (symsize) eq 0) then symsize= charsize + intarr(n)
+if n_elements(number) eq 0 then number = 1
+ case N_elements(colorsi) of
+ 0: colors = replicate(!P.color,n) ;Default is !P.COLOR
+ 1: colors = replicate(colorsi,n)
+ else: colors = colorsi
+ endcase
+
+ case N_elements(textcolorsi) of
+ 0: textcolors = replicate(!P.color,n) ;Default is !P.COLOR
+ 1: textcolors = replicate(textcolorsi,n)
+ else: textcolors = textcolorsi
+ endcase
+ fill = keyword_set(fill)
+if n_elements(usersym) eq 1 then usersym = 2*[[0,0],[0,1],[1,1],[1,0],[0,0]]-1
+;
+; =====>> INITIALIZE SPACING
+;
+if n_elements(spacing) eq 0 then spacing = 1.2
+if n_elements(pspacing) eq 0 then pspacing = 3
+xspacing = !d.x_ch_size/float(!d.x_size) * (spacing > charsize)
+yspacing = !d.y_ch_size/float(!d.y_size) * (spacing > charsize)
+ltor = 1 ; flag for left-to-right
+if n_elements(left) eq 1 then ltor = left eq 1
+if n_elements(right) eq 1 then ltor = right ne 1
+ttob = 1 ; flag for top-to-bottom
+if n_elements(top) eq 1 then ttob = top eq 1
+if n_elements(bottom) eq 1 then ttob = bottom ne 1
+xalign = ltor ne 1 ; x alignment: 1 or 0
+yalign = -0.5*ttob + 1 ; y alignment: 0.5 or 1
+xsign = 2*ltor - 1 ; xspacing direction: 1 or -1
+ysign = 2*ttob - 1 ; yspacing direction: 1 or -1
+if not ttob then yspacing = -yspacing
+if not ltor then xspacing = -xspacing
+;
+; =====>> INITIALIZE POSITIONS: FIRST CALCULATE X OFFSET FOR TEXT
+;
+xt = 0
+if nlpv gt 0 then begin ; SKIP IF TEXT ITEMS ONLY.
+if vertical then begin ; CALC OFFSET FOR TEXT START
+ for i = 0,n-1 do begin
+ if (psym[i] eq 0) and (vectorfont[i] eq '') then num = (number + 1) > 3 else num = number
+ if psym[i] lt 0 then num = number > 2 ; TO SHOW CONNECTING LINE
+ if psym[i] eq 0 then expand = 1 else expand = 2
+ thisxt = (expand*pspacing*(num-1)*xspacing)
+ if ltor then xt = thisxt > xt else xt = thisxt < xt
+ endfor
+endif ; NOW xt IS AN X OFFSET TO ALIGN ALL TEXT ENTRIES.
+endif
+;
+; =====>> INITIALIZE POSITIONS: SECOND LOCATE BORDER
+;
+if !x.window[0] eq !x.window[1] then begin
+ plot,/nodata,xstyle=4,ystyle=4,[0],/noerase
+endif
+; next line takes care of weirdness with small windows
+pos = [min(!x.window),min(!y.window),max(!x.window),max(!y.window)]
+case n_elements(position) of
+ 0: begin
+ if ltor then px = pos[0] else px = pos[2]
+ if ttob then py = pos[3] else py = pos[1]
+ if keyword_set(center) then begin
+ if not keyword_set(right) and not keyword_set(left) then $
+ px = (pos[0] + pos[2])/2. - xt
+ if not keyword_set(top) and not keyword_set(bottom) then $
+ py = (pos[1] + pos[3])/2. + n*yspacing
+ endif
+ position = [px,py] + [xspacing,-yspacing]
+ end
+ 1: begin ; interactive
+ message,/inform,'Place mouse at upper left corner and click any mouse button.'
+ cursor,x,y,/normal
+ position = [x,y]
+ end
+ 2: begin ; convert upper left corner to normal coordinates
+ if keyword_set(data) then $
+ position = convert_coord(position,/to_norm) $
+ else if keyword_set(device) then $
+ position = convert_coord(position,/to_norm,/device) $
+ else if not keyword_set(normal) then $
+ position = convert_coord(position,/to_norm)
+ end
+ else: message,'Position keyword can have 0, 1, or 2 elements only. Try legend,/help.'
+endcase
+
+yoff = 0.25*yspacing*ysign ; VERT. OFFSET FOR SYM/LINE.
+
+x0 = position[0] + (margin)*xspacing ; INITIAL X & Y POSITIONS
+y0 = position[1] - margin*yspacing + yalign*yspacing ; WELL, THIS WORKS!
+;
+; =====>> OUTPUT TEXT FOR LEGEND, ITEM BY ITEM.
+; =====>> FOR EACH ITEM, PLACE SYM/LINE, THEN DELIMITER,
+; =====>> THEN TEXT---UPDATING X & Y POSITIONS EACH TIME.
+; =====>> THERE ARE A NUMBER OF EXCEPTIONS DONE WITH IF STATEMENTS.
+;
+for iclr = 0,clear do begin
+ y = y0 ; STARTING X & Y POSITIONS
+ x = x0
+ if ltor then xend = 0 else xend = 1 ; SAVED WIDTH FOR DRAWING BOX
+
+ if ttob then ii = [0,n-1,1] else ii = [n-1,0,-1]
+ for i = ii[0],ii[1],ii[2] do begin
+ if vertical then x = x0 else y = y0 ; RESET EITHER X OR Y
+ x = x + xspacing ; UPDATE X & Y POSITIONS
+ y = y - yspacing
+ if nlpv eq 0 then goto,TEXT_ONLY ; FLAG FOR TEXT ONLY
+ if (psym[i] eq 0) and (vectorfont[i] eq '') then num = (number + 1) > 3 else num = number
+ if psym[i] lt 0 then num = number > 2 ; TO SHOW CONNECTING LINE
+ if psym[i] eq 0 then expand = 1 else expand = 2
+ xp = x + expand*pspacing*indgen(num)*xspacing
+ if (psym[i] gt 0) and (num eq 1) and vertical then xp = x + xt/2.
+ yp = y + intarr(num)
+ if vectorfont[i] eq '' then yp = yp + yoff
+ if psym[i] eq 0 then begin
+ xp = [min(xp),max(xp)] ; TO EXPOSE LINESTYLES
+ yp = [min(yp),max(yp)] ; DITTO
+ endif
+ if (psym[i] eq 8) and (N_elements(usersym) GT 1) then $
+ usersym,usersym,fill=fill,color=colors[i]
+;; extra by djseed .. psym=88 means use the already defined usersymbol
+ if psym[i] eq 88 then psym[i] =8
+ if vectorfont[i] ne '' then begin
+; if (num eq 1) and vertical then xp = x + xt/2 ; IF 1, CENTERED.
+ xyouts,xp,yp,vectorfont[i],width=width,color=colors[i] $
+ ,size=charsize,align=xalign,charthick = charthick,/norm
+ xt = xt > width
+ xp = xp + width/2.
+ endif else begin
+ if symline and (linestyle[i] ge 0) then plots,xp,yp,color=colors[i] $
+ ,/normal,linestyle=linestyle[i],psym=psym[i],symsize=symsize[i], $
+ thick=thick[i]
+ endelse
+
+ if vertical then x = x + xt else if ltor then x = max(xp) else x = min(xp)
+ if symline then x = x + xspacing
+ TEXT_ONLY:
+ if vertical and (vectorfont[i] eq '') and symline and (linestyle[i] eq -99) then x=x0 + xspacing
+ xyouts,x,y,delimiter,width=width,/norm,color=textcolors[i], $
+ size=charsize,align=xalign,charthick = charthick
+ x = x + width*xsign
+ if width ne 0 then x = x + 0.5*xspacing
+ xyouts,x,y,items[i],width=width,/norm,color=textcolors[i],size=charsize, $
+ align=xalign,charthick=charthick
+ x = x + width*xsign
+ if not vertical and (i lt (n-1)) then x = x+2*xspacing; ADD INTER-ITEM SPACE
+ xfinal = (x + xspacing*margin)
+ if ltor then xend = xfinal > xend else xend = xfinal < xend ; UPDATE END X
+ endfor
+
+ if (iclr lt clear ) then begin
+; =====>> CLEAR AREA
+ x = position[0]
+ y = position[1]
+ if vertical then bottom = n else bottom = 1
+ ywidth = - (2*margin+bottom-0.5)*yspacing
+ corners = [x,y+ywidth,xend,y]
+ polyfill,[x,xend,xend,x,x],y + [0,0,ywidth,ywidth,0],/norm,color=-1
+; plots,[x,xend,xend,x,x],y + [0,0,ywidth,ywidth,0],thick=2
+ endif else begin
+
+;
+; =====>> OUTPUT BORDER
+;
+ x = position[0]
+ y = position[1]
+ if vertical then bottom = n else bottom = 1
+ ywidth = - (2*margin+bottom-0.5)*yspacing
+ corners = [x,y+ywidth,xend,y]
+ if box then plots,[x,xend,xend,x,x],y + [0,0,ywidth,ywidth,0],/norm
+ return
+ endelse
+endfor
+
+end
+
diff --git a/examples/Disk-Patch/GravityOnly/makeIC.py b/examples/Disk-Patch/GravityOnly/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..135d1d4e34e0ec4c32d75a2c49d6e55f46db4e1c
--- /dev/null
+++ b/examples/Disk-Patch/GravityOnly/makeIC.py
@@ -0,0 +1,161 @@
+###############################################################################
+ # 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 box of [0:BoxSize,0:BoxSize,-2scale_height:2scale_height]
+# see Creasey, Theuns & Bower, 2013, for the equations:
+# disk parameters are: surface density sigma
+# scale height b
+# density: rho(z) = (sigma/2b) sech^2(z/b)
+# isothermal velocity dispersion = <v_z^2? = b pi G sigma
+# grad potential = 2 pi G sigma tanh(z/b)
+# potential = ln(cosh(z/b)) + const
+# Dynamical time = sqrt(b / (G sigma))
+# to obtain the 1/ch^2(z/b) profile from a uniform profile (a glass, say, or a uniform random variable), note that, when integrating in z
+# \int 0^z dz/ch^2(z) = tanh(z)-tanh(0) = \int_0^x dx = x (where the last integral refers to a uniform density distribution), so that z = atanh(x)
+# usage: python makeIC.py 1000
+
+# 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 = (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
+
+
+# parameters of potential
+surface_density = 10.
+scale_height = 100.
+
+# 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
+v_disp = numpy.sqrt(scale_height * math.pi * const_G * surface_density)
+t_dyn = numpy.sqrt(scale_height / (const_G * surface_density))
+print 'dynamical time = ',t_dyn
+print ' velocity dispersion = ',v_disp
+
+# Parameters
+periodic= 1 # 1 For periodic box
+boxSize = 600. #
+Radius = 100. # maximum radius of particles [kpc]
+G = const_G
+
+N = int(sys.argv[1]) # Number of particles
+
+# these are not used but necessary for I/O
+rho = 2. # Density
+P = 1. # Pressure
+gamma = 5./3. # Gas adiabatic index
+fileName = "Disk-Patch.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
+r = numpy.zeros((numPart, 3))
+r[:,0] = numpy.random.rand(N) * boxSize
+r[:,1] = numpy.random.rand(N) * boxSize
+z = scale_height * numpy.arctanh(numpy.random.rand(2*N))
+gd = z < boxSize / 2
+r[:,2] = z[gd][0:N]
+random = numpy.random.rand(N) > 0.5
+r[random,2] *= -1
+r[:,2] += 0.5 * boxSize
+
+#generate particle velocities
+v = numpy.zeros((numPart, 3))
+v = numpy.zeros(1)
+#v[:,2] =
+
+
+ds = grp1.create_dataset('Velocities', (numPart, 3), 'f')
+ds[()] = v
+
+
+m = numpy.ones((numPart, ), dtype=numpy.float32) * mass
+ds = grp1.create_dataset('Masses', (numPart,), 'f')
+ds[()] = m
+m = 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()
diff --git a/examples/Disk-Patch/GravityOnly/readme.txt b/examples/Disk-Patch/GravityOnly/readme.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a18365146aa13abb81e7d3b8fcac31b92959ae3f
--- /dev/null
+++ b/examples/Disk-Patch/GravityOnly/readme.txt
@@ -0,0 +1,13 @@
+sets-up for a potential of a patch f disk, see Creasey, Theuns &
+Bower,2013.
+
+rho(z) = (Sigma/2b) / cosh^2(z/b)
+where Sigma is the surface density, and b the scale height
+the corresponding force is
+dphi/dz = 2 pi G Sigma tanh(z/b)
+
+which satifies d^2phi/dz^2 = 4 pi G rho
+
+
+
+
diff --git a/examples/Disk-Patch/GravityOnly/run.sh b/examples/Disk-Patch/GravityOnly/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..89e7c10ce19ce0576127551f63eea5c9e4e955ff
--- /dev/null
+++ b/examples/Disk-Patch/GravityOnly/run.sh
@@ -0,0 +1,10 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e Isothermal.hdf5 ]
+then
+ echo "Generating initial conditions for the point mass potential box example..."
+ python makeIC.py 1000
+fi
+
+../swift -g -t 2 externalGravity.yml
diff --git a/examples/Disk-Patch/GravityOnly/test.pro b/examples/Disk-Patch/GravityOnly/test.pro
new file mode 100644
index 0000000000000000000000000000000000000000..4bd0d001975d80b6729cf2ef7b95f81da5bc4fe8
--- /dev/null
+++ b/examples/Disk-Patch/GravityOnly/test.pro
@@ -0,0 +1,158 @@
+;
+; 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 = 'Disk-Patch_'
+
+; set properties of potential
+uL = phys.pc ; unit of length
+uM = phys.msun ; unit of mass
+uV = 1d5 ; unit of velocity
+
+; properties of patch
+surface_density = 10.
+scale_height = 100.
+
+; derived units
+constG = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ;
+pcentre = [0.,0.,300.] * 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) ; 2 pi G Sigma b ln(cosh(z/b)) + const
+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)
+
+;; ; if you want to sort particles by 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,*]
+ dz = reform(p[2,0:nfollow-1]-pcentre[2])
+; print,'time = ',time,p[0,0],v[0,0],id[0]
+ ek = 0.5 * dv
+ ep = fltarr(nfollow)
+ ep = 2 * !pi * constG * surface_density * scale_height * alog(cosh(abs(dz)/scale_height))
+ 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]
+ 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 vertical oscillation
+ win,2
+ xr = [min(tout), max(tout)]
+ yr = [-3,3]*scale_height
+ 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,tout,zout[i,*],color=color(i)
+ screen_to_png,'orbit.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
+
+
diff --git a/examples/Disk-Patch/HydroStatic/disk-patch.yml b/examples/Disk-Patch/HydroStatic/disk-patch.yml
new file mode 100644
index 0000000000000000000000000000000000000000..46cef05130eb337de25d15d9f34e2c346a9e314e
--- /dev/null
+++ b/examples/Disk-Patch/HydroStatic/disk-patch.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: 480. # The end time of the simulation (in internal units).
+ dt_min: 1e-5 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 4. # 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: 8. # 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/legend.pro b/examples/Disk-Patch/HydroStatic/legend.pro
new file mode 100755
index 0000000000000000000000000000000000000000..e510bdfcd7dcabed796cb090c6fa7d54536608b6
--- /dev/null
+++ b/examples/Disk-Patch/HydroStatic/legend.pro
@@ -0,0 +1,459 @@
+;+
+; NAME:
+; LEGEND
+; PURPOSE:
+; Create an annotation legend for a plot.
+; EXPLANATION:
+; This procedure makes a legend for a plot. The legend can contain
+; a mixture of symbols, linestyles, Hershey characters (vectorfont),
+; and filled polygons (usersym). A test procedure, legendtest.pro,
+; shows legend's capabilities. Placement of the legend is controlled
+; with keywords like /right, /top, and /center or by using a position
+; keyword for exact placement (position=[x,y]) or via mouse (/position).
+; CALLING SEQUENCE:
+; LEGEND [,items][,keyword options]
+; EXAMPLES:
+; The call:
+; legend,['Plus sign','Asterisk','Period'],psym=[1,2,3]
+; produces:
+; -----------------
+; | |
+; | + Plus sign |
+; | * Asterisk |
+; | . Period |
+; | |
+; -----------------
+; Each symbol is drawn with a plots command, so they look OK.
+; Other examples are given in optional output keywords.
+;
+; lines = indgen(6) ; for line styles
+; items = 'linestyle '+strtrim(lines,2) ; annotations
+; legend,items,linestyle=lines ; vertical legend---upper left
+; items = ['Plus sign','Asterisk','Period']
+; sym = [1,2,3]
+; legend,items,psym=sym ; ditto except using symbols
+; legend,items,psym=sym,/horizontal ; horizontal format
+; legend,items,psym=sym,box=0 ; sans border
+; legend,items,psym=sym,delimiter='=' ; embed '=' betw psym & text
+; legend,items,psym=sym,margin=2 ; 2-character margin
+; legend,items,psym=sym,position=[x,y] ; upper left in data coords
+; legend,items,psym=sym,pos=[x,y],/norm ; upper left in normal coords
+; legend,items,psym=sym,pos=[x,y],/device ; upper left in device coords
+; legend,items,psym=sym,/position ; interactive position
+; legend,items,psym=sym,/right ; at upper right
+; legend,items,psym=sym,/bottom ; at lower left
+; legend,items,psym=sym,/center ; approximately near center
+; legend,items,psym=sym,number=2 ; plot two symbols, not one
+; legend,items,/fill,psym=[8,8,8],colors=[10,20,30]; 3 filled squares
+; INPUTS:
+; items = text for the items in the legend, a string array.
+; For example, items = ['diamond','asterisk','square'].
+; You can omit items if you don't want any text labels.
+; OPTIONAL INPUT KEYWORDS:
+;
+; linestyle = array of linestyle numbers If linestyle[i] < 0, then omit
+; ith symbol or line to allow a multi-line entry. If
+; linestyle = -99 then text will be left-justified.
+; psym = array of plot symbol numbers. If psym[i] is negative, then a
+; line connects pts for ith item. If psym[i] = 8, then the
+; procedure usersym is called with vertices define in the
+; keyword usersym. If psym[i] = 88, then use the previously
+; defined user symbol
+; vectorfont = vector-drawn characters for the sym/line column, e.g.,
+; ['!9B!3','!9C!3','!9D!3'] produces an open square, a checkmark,
+; and a partial derivative, which might have accompanying items
+; ['BOX','CHECK','PARTIAL DERIVATIVE'].
+; There is no check that !p.font is set properly, e.g., -1 for
+; X and 0 for PostScript. This can produce an error, e.g., use
+; !20 with PostScript and !p.font=0, but allows use of Hershey
+; *AND* PostScript fonts together.
+; N. B.: Choose any of linestyle, psym, and/or vectorfont. If none is
+; present, only the text is output. If more than one
+; is present, all need the same number of elements, and normal
+; plot behaviour occurs.
+; By default, if psym is positive, you get one point so there is
+; no connecting line. If vectorfont[i] = '',
+; then plots is called to make a symbol or a line, but if
+; vectorfont[i] is a non-null string, then xyouts is called.
+; /help = flag to print header
+; /horizontal = flag to make the legend horizontal
+; /vertical = flag to make the legend vertical (D=vertical)
+; box = flag to include/omit box around the legend (D=include)
+; clear = flag to clear the box area before drawing the legend
+; delimiter = embedded character(s) between symbol and text (D=none)
+; colors = array of colors for plot symbols/lines (D=!P.color)
+; textcolors = array of colors for text (D=!P.color)
+; margin = margin around text measured in characters and lines
+; spacing = line spacing (D=bit more than character height)
+; pspacing = psym spacing (D=3 characters) (when number of symbols is
+; greater than 1)
+; charsize = just like !p.charsize for plot labels
+; charthick = just like !p.charthick for plot labels
+; thick = array of line thickness numbers (D = !P.thick), if used, then
+; linestyle must also be specified
+; position = data coordinates of the /top (D) /left (D) of the legend
+; normal = use normal coordinates for position, not data
+; device = use device coordinates for position, not data
+; number = number of plot symbols to plot or length of line (D=1)
+; usersym = 2-D array of vertices, cf. usersym in IDL manual.
+; (/USERSYM =square, default is to use existing USERSYM definition)
+; /fill = flag to fill the usersym
+; /left_legend = flag to place legend snug against left side of plot
+; window (D)
+; /right_legend = flag to place legend snug against right side of plot
+; window. If /right,pos=[x,y], then x is position of RHS and
+; text runs right-to-left.
+; /top_legend = flag to place legend snug against top of plot window (D)
+; /bottom = flag to place legend snug against bottom of plot window
+; /top,pos=[x,y] and /bottom,pos=[x,y] produce same positions.
+;
+; If LINESTYLE, PSYM, VECTORFONT, THICK, COLORS, or TEXTCOLORS are
+; supplied as scalars, then the scalar value is set for every line or
+; symbol in the legend.
+; Outputs:
+; legend to current plot device
+; OPTIONAL OUTPUT KEYWORDS:
+; corners = 4-element array, like !p.position, of the normalized
+; coords for the box (even if box=0): [llx,lly,urx,ury].
+; Useful for multi-column or multi-line legends, for example,
+; to make a 2-column legend, you might do the following:
+; c1_items = ['diamond','asterisk','square']
+; c1_psym = [4,2,6]
+; c2_items = ['solid','dashed','dotted']
+; c2_line = [0,2,1]
+; legend,c1_items,psym=c1_psym,corners=c1,box=0
+; legend,c2_items,line=c2_line,corners=c2,box=0,pos=[c1[2],c1[3]]
+; c = [c1[0]<c2[0],c1[1]<c2[1],c1[2]>c2[2],c1[3]>c2[3]]
+; plots,[c[0],c[0],c[2],c[2],c[0]],[c[1],c[3],c[3],c[1],c[1]],/norm
+; Useful also to place the legend. Here's an automatic way to place
+; the legend in the lower right corner. The difficulty is that the
+; legend's width is unknown until it is plotted. In this example,
+; the legend is plotted twice: the first time in the upper left, the
+; second time in the lower right.
+; legend,['1','22','333','4444'],linestyle=indgen(4),corners=corners
+; ; BOGUS LEGEND---FIRST TIME TO REPORT CORNERS
+; xydims = [corners[2]-corners[0],corners[3]-corners[1]]
+; ; SAVE WIDTH AND HEIGHT
+; chdim=[!d.x_ch_size/float(!d.x_size),!d.y_ch_size/float(!d.y_size)]
+; ; DIMENSIONS OF ONE CHARACTER IN NORMALIZED COORDS
+; pos = [!x.window[1]-chdim[0]-xydims[0] $
+; ,!y.window[0]+chdim[1]+xydims[1]]
+; ; CALCULATE POSITION FOR LOWER RIGHT
+; plot,findgen(10) ; SIMPLE PLOT; YOU DO WHATEVER YOU WANT HERE.
+; legend,['1','22','333','4444'],linestyle=indgen(4),pos=pos
+; ; REDO THE LEGEND IN LOWER RIGHT CORNER
+; You can modify the pos calculation to place the legend where you
+; want. For example to place it in the upper right:
+; pos = [!x.window[1]-chdim[0]-xydims[0],!y.window[1]-xydims[1]]
+; Common blocks:
+; none
+; Procedure:
+; If keyword help is set, call doc_library to print header.
+; See notes in the code. Much of the code deals with placement of the
+; legend. The main problem with placement is not being
+; able to sense the length of a string before it is output. Some crude
+; approximations are used for centering.
+; Restrictions:
+; Here are some things that aren't implemented.
+; - An orientation keyword would allow lines at angles in the legend.
+; - An array of usersyms would be nice---simple change.
+; - An order option to interchange symbols and text might be nice.
+; - Somebody might like double boxes, e.g., with box = 2.
+; - Another feature might be a continuous bar with ticks and text.
+; - There are no guards to avoid writing outside the plot area.
+; - There is no provision for multi-line text, e.g., '1st line!c2nd line'
+; Sensing !c would be easy, but !c isn't implemented for PostScript.
+; A better way might be to simply output the 2nd line as another item
+; but without any accompanying symbol or linestyle. A flag to omit
+; the symbol and linestyle is linestyle[i] = -1.
+; - There is no ability to make a title line containing any of titles
+; for the legend, for the symbols, or for the text.
+; Side Effects:
+; Modification history:
+; write, 24-25 Aug 92, F K Knight (knight@ll.mit.edu)
+; allow omission of items or omission of both psym and linestyle, add
+; corners keyword to facilitate multi-column legends, improve place-
+; ment of symbols and text, add guards for unequal size, 26 Aug 92, FKK
+; add linestyle(i)=-1 to suppress a single symbol/line, 27 Aug 92, FKK
+; add keyword vectorfont to allow characters in the sym/line column,
+; 28 Aug 92, FKK
+; add /top, /bottom, /left, /right keywords for automatic placement at
+; the four corners of the plot window. The /right keyword forces
+; right-to-left printing of menu. 18 Jun 93, FKK
+; change default position to data coords and add normal, data, and
+; device keywords, 17 Jan 94, FKK
+; add /center keyword for positioning, but it is not precise because
+; text string lengths cannot be known in advance, 17 Jan 94, FKK
+; add interactive positioning with /position keyword, 17 Jan 94, FKK
+; allow a legend with just text, no plotting symbols. This helps in
+; simply describing a plot or writing assumptions done, 4 Feb 94, FKK
+; added thick, symsize, and clear keyword Feb 96, W. Landsman HSTX
+; David Seed, HR Wallingford, d.seed@hrwallingford.co.uk
+; allow scalar specification of keywords, Mar 96, W. Landsman HSTX
+; added charthick keyword, June 96, W. Landsman HSTX
+; Made keyword names left,right,top,bottom,center longer,
+; Aug 16, 2000, Kim Tolbert
+; Added ability to have regular text lines in addition to plot legend
+; lines in legend. If linestyle is -99 that item is left-justified.
+; Previously, only option for no sym/line was linestyle=-1, but then text
+; was lined up after sym/line column. 10 Oct 2000, Kim Tolbert
+; Make default value of thick = !P.thick W. Landsman Jan. 2001
+; Don't overwrite existing USERSYM definition W. Landsman Mar. 2002
+;-
+pro legend, items, BOTTOM_LEGEND=bottom, BOX = box, CENTER_LEGEND=center, $
+ CHARTHICK=charthick, CHARSIZE = charsize, CLEAR = clear, COLORS = colorsi, $
+ CORNERS = corners, DATA=data, DELIMITER=delimiter, DEVICE=device, $
+ FILL=fill, HELP = help, HORIZONTAL=horizontal,LEFT_LEGEND=left, $
+ LINESTYLE=linestylei, MARGIN=margin, NORMAL=normal, NUMBER=number, $
+ POSITION=position,PSPACING=pspacing, PSYM=psymi, RIGHT_LEGEND=right, $
+ SPACING=spacing, SYMSIZE=symsize, TEXTCOLORS=textcolorsi, THICK=thicki, $
+ TOP_LEGEND=top, USERSYM=usersym, VECTORFONT=vectorfonti, VERTICAL=vertical
+;
+; =====>> HELP
+;
+on_error,2
+if keyword_set(help) then begin & doc_library,'legend' & return & endif
+;
+; =====>> SET DEFAULTS FOR SYMBOLS, LINESTYLES, AND ITEMS.
+;
+ ni = n_elements(items)
+ np = n_elements(psymi)
+ nl = n_elements(linestylei)
+ nth = n_elements(thicki)
+ nv = n_elements(vectorfonti)
+ nlpv = max([np,nl,nv])
+ n = max([ni,np,nl,nv]) ; NUMBER OF ENTRIES
+strn = strtrim(n,2) ; FOR ERROR MESSAGES
+if n eq 0 then message,'No inputs! For help, type legend,/help.'
+if ni eq 0 then begin
+ items = replicate('',n) ; DEFAULT BLANK ARRAY
+endif else begin
+ if size(items,/TNAME) NE 'STRING' then message, $
+ 'First parameter must be a string array. For help, type legend,/help.'
+ if ni ne n then message,'Must have number of items equal to '+strn
+endelse
+symline = (np ne 0) or (nl ne 0) ; FLAG TO PLOT SYM/LINE
+ if (np ne 0) and (np ne n) and (np NE 1) then message, $
+ 'Must have 0, 1 or '+strn+' elements in PSYM array.'
+ if (nl ne 0) and (nl ne n) and (nl NE 1) then message, $
+ 'Must have 0, 1 or '+strn+' elements in LINESTYLE array.'
+ if (nth ne 0) and (nth ne n) and (nth NE 1) then message, $
+ 'Must have 0, 1 or '+strn+' elements in THICK array.'
+
+ case nl of
+ 0: linestyle = intarr(n) ;Default = solid
+ 1: linestyle = intarr(n) + linestylei
+ else: linestyle = linestylei
+ endcase
+
+ case nth of
+ 0: thick = replicate(!p.thick,n) ;Default = !P.THICK
+ 1: thick = intarr(n) + thicki
+ else: thick = thicki
+ endcase
+
+ case np of ;Get symbols
+ 0: psym = intarr(n) ;Default = solid
+ 1: psym = intarr(n) + psymi
+ else: psym = psymi
+ endcase
+
+ case nv of
+ 0: vectorfont = replicate('',n)
+ 1: vectorfont = replicate(vectorfonti,n)
+ else: vectorfont = vectorfonti
+ endcase
+;
+; =====>> CHOOSE VERTICAL OR HORIZONTAL ORIENTATION.
+;
+if n_elements(horizontal) eq 0 then begin ; D=VERTICAL
+ if n_elements(vertical) eq 0 then vertical = 1
+endif else begin
+ if n_elements(vertical) eq 0 then vertical = not horizontal
+endelse
+;
+; =====>> SET DEFAULTS FOR OTHER OPTIONS.
+;
+if n_elements(box) eq 0 then box = 1
+if n_elements(clear) eq 0 then clear = 0
+
+if n_elements(margin) eq 0 then margin = 0.5
+if n_elements(delimiter) eq 0 then delimiter = ''
+if n_elements(charsize) eq 0 then charsize = !p.charsize
+if n_elements(charthick) eq 0 then charthick = !p.charthick
+if charsize eq 0 then charsize = 1
+if (n_elements (symsize) eq 0) then symsize= charsize + intarr(n)
+if n_elements(number) eq 0 then number = 1
+ case N_elements(colorsi) of
+ 0: colors = replicate(!P.color,n) ;Default is !P.COLOR
+ 1: colors = replicate(colorsi,n)
+ else: colors = colorsi
+ endcase
+
+ case N_elements(textcolorsi) of
+ 0: textcolors = replicate(!P.color,n) ;Default is !P.COLOR
+ 1: textcolors = replicate(textcolorsi,n)
+ else: textcolors = textcolorsi
+ endcase
+ fill = keyword_set(fill)
+if n_elements(usersym) eq 1 then usersym = 2*[[0,0],[0,1],[1,1],[1,0],[0,0]]-1
+;
+; =====>> INITIALIZE SPACING
+;
+if n_elements(spacing) eq 0 then spacing = 1.2
+if n_elements(pspacing) eq 0 then pspacing = 3
+xspacing = !d.x_ch_size/float(!d.x_size) * (spacing > charsize)
+yspacing = !d.y_ch_size/float(!d.y_size) * (spacing > charsize)
+ltor = 1 ; flag for left-to-right
+if n_elements(left) eq 1 then ltor = left eq 1
+if n_elements(right) eq 1 then ltor = right ne 1
+ttob = 1 ; flag for top-to-bottom
+if n_elements(top) eq 1 then ttob = top eq 1
+if n_elements(bottom) eq 1 then ttob = bottom ne 1
+xalign = ltor ne 1 ; x alignment: 1 or 0
+yalign = -0.5*ttob + 1 ; y alignment: 0.5 or 1
+xsign = 2*ltor - 1 ; xspacing direction: 1 or -1
+ysign = 2*ttob - 1 ; yspacing direction: 1 or -1
+if not ttob then yspacing = -yspacing
+if not ltor then xspacing = -xspacing
+;
+; =====>> INITIALIZE POSITIONS: FIRST CALCULATE X OFFSET FOR TEXT
+;
+xt = 0
+if nlpv gt 0 then begin ; SKIP IF TEXT ITEMS ONLY.
+if vertical then begin ; CALC OFFSET FOR TEXT START
+ for i = 0,n-1 do begin
+ if (psym[i] eq 0) and (vectorfont[i] eq '') then num = (number + 1) > 3 else num = number
+ if psym[i] lt 0 then num = number > 2 ; TO SHOW CONNECTING LINE
+ if psym[i] eq 0 then expand = 1 else expand = 2
+ thisxt = (expand*pspacing*(num-1)*xspacing)
+ if ltor then xt = thisxt > xt else xt = thisxt < xt
+ endfor
+endif ; NOW xt IS AN X OFFSET TO ALIGN ALL TEXT ENTRIES.
+endif
+;
+; =====>> INITIALIZE POSITIONS: SECOND LOCATE BORDER
+;
+if !x.window[0] eq !x.window[1] then begin
+ plot,/nodata,xstyle=4,ystyle=4,[0],/noerase
+endif
+; next line takes care of weirdness with small windows
+pos = [min(!x.window),min(!y.window),max(!x.window),max(!y.window)]
+case n_elements(position) of
+ 0: begin
+ if ltor then px = pos[0] else px = pos[2]
+ if ttob then py = pos[3] else py = pos[1]
+ if keyword_set(center) then begin
+ if not keyword_set(right) and not keyword_set(left) then $
+ px = (pos[0] + pos[2])/2. - xt
+ if not keyword_set(top) and not keyword_set(bottom) then $
+ py = (pos[1] + pos[3])/2. + n*yspacing
+ endif
+ position = [px,py] + [xspacing,-yspacing]
+ end
+ 1: begin ; interactive
+ message,/inform,'Place mouse at upper left corner and click any mouse button.'
+ cursor,x,y,/normal
+ position = [x,y]
+ end
+ 2: begin ; convert upper left corner to normal coordinates
+ if keyword_set(data) then $
+ position = convert_coord(position,/to_norm) $
+ else if keyword_set(device) then $
+ position = convert_coord(position,/to_norm,/device) $
+ else if not keyword_set(normal) then $
+ position = convert_coord(position,/to_norm)
+ end
+ else: message,'Position keyword can have 0, 1, or 2 elements only. Try legend,/help.'
+endcase
+
+yoff = 0.25*yspacing*ysign ; VERT. OFFSET FOR SYM/LINE.
+
+x0 = position[0] + (margin)*xspacing ; INITIAL X & Y POSITIONS
+y0 = position[1] - margin*yspacing + yalign*yspacing ; WELL, THIS WORKS!
+;
+; =====>> OUTPUT TEXT FOR LEGEND, ITEM BY ITEM.
+; =====>> FOR EACH ITEM, PLACE SYM/LINE, THEN DELIMITER,
+; =====>> THEN TEXT---UPDATING X & Y POSITIONS EACH TIME.
+; =====>> THERE ARE A NUMBER OF EXCEPTIONS DONE WITH IF STATEMENTS.
+;
+for iclr = 0,clear do begin
+ y = y0 ; STARTING X & Y POSITIONS
+ x = x0
+ if ltor then xend = 0 else xend = 1 ; SAVED WIDTH FOR DRAWING BOX
+
+ if ttob then ii = [0,n-1,1] else ii = [n-1,0,-1]
+ for i = ii[0],ii[1],ii[2] do begin
+ if vertical then x = x0 else y = y0 ; RESET EITHER X OR Y
+ x = x + xspacing ; UPDATE X & Y POSITIONS
+ y = y - yspacing
+ if nlpv eq 0 then goto,TEXT_ONLY ; FLAG FOR TEXT ONLY
+ if (psym[i] eq 0) and (vectorfont[i] eq '') then num = (number + 1) > 3 else num = number
+ if psym[i] lt 0 then num = number > 2 ; TO SHOW CONNECTING LINE
+ if psym[i] eq 0 then expand = 1 else expand = 2
+ xp = x + expand*pspacing*indgen(num)*xspacing
+ if (psym[i] gt 0) and (num eq 1) and vertical then xp = x + xt/2.
+ yp = y + intarr(num)
+ if vectorfont[i] eq '' then yp = yp + yoff
+ if psym[i] eq 0 then begin
+ xp = [min(xp),max(xp)] ; TO EXPOSE LINESTYLES
+ yp = [min(yp),max(yp)] ; DITTO
+ endif
+ if (psym[i] eq 8) and (N_elements(usersym) GT 1) then $
+ usersym,usersym,fill=fill,color=colors[i]
+;; extra by djseed .. psym=88 means use the already defined usersymbol
+ if psym[i] eq 88 then psym[i] =8
+ if vectorfont[i] ne '' then begin
+; if (num eq 1) and vertical then xp = x + xt/2 ; IF 1, CENTERED.
+ xyouts,xp,yp,vectorfont[i],width=width,color=colors[i] $
+ ,size=charsize,align=xalign,charthick = charthick,/norm
+ xt = xt > width
+ xp = xp + width/2.
+ endif else begin
+ if symline and (linestyle[i] ge 0) then plots,xp,yp,color=colors[i] $
+ ,/normal,linestyle=linestyle[i],psym=psym[i],symsize=symsize[i], $
+ thick=thick[i]
+ endelse
+
+ if vertical then x = x + xt else if ltor then x = max(xp) else x = min(xp)
+ if symline then x = x + xspacing
+ TEXT_ONLY:
+ if vertical and (vectorfont[i] eq '') and symline and (linestyle[i] eq -99) then x=x0 + xspacing
+ xyouts,x,y,delimiter,width=width,/norm,color=textcolors[i], $
+ size=charsize,align=xalign,charthick = charthick
+ x = x + width*xsign
+ if width ne 0 then x = x + 0.5*xspacing
+ xyouts,x,y,items[i],width=width,/norm,color=textcolors[i],size=charsize, $
+ align=xalign,charthick=charthick
+ x = x + width*xsign
+ if not vertical and (i lt (n-1)) then x = x+2*xspacing; ADD INTER-ITEM SPACE
+ xfinal = (x + xspacing*margin)
+ if ltor then xend = xfinal > xend else xend = xfinal < xend ; UPDATE END X
+ endfor
+
+ if (iclr lt clear ) then begin
+; =====>> CLEAR AREA
+ x = position[0]
+ y = position[1]
+ if vertical then bottom = n else bottom = 1
+ ywidth = - (2*margin+bottom-0.5)*yspacing
+ corners = [x,y+ywidth,xend,y]
+ polyfill,[x,xend,xend,x,x],y + [0,0,ywidth,ywidth,0],/norm,color=-1
+; plots,[x,xend,xend,x,x],y + [0,0,ywidth,ywidth,0],thick=2
+ endif else begin
+
+;
+; =====>> OUTPUT BORDER
+;
+ x = position[0]
+ y = position[1]
+ if vertical then bottom = n else bottom = 1
+ ywidth = - (2*margin+bottom-0.5)*yspacing
+ corners = [x,y+ywidth,xend,y]
+ if box then plots,[x,xend,xend,x,x],y + [0,0,ywidth,ywidth,0],/norm
+ return
+ endelse
+endfor
+
+end
+
diff --git a/examples/Disk-Patch/HydroStatic/makeIC-stretch.py b/examples/Disk-Patch/HydroStatic/makeIC-stretch.py
new file mode 100644
index 0000000000000000000000000000000000000000..761f50eec157da6a1841e885fb89944f8fc41ba6
--- /dev/null
+++ b/examples/Disk-Patch/HydroStatic/makeIC-stretch.py
@@ -0,0 +1,285 @@
+###############################################################################
+ # 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
+import matplotlib.pyplot as plt
+
+# Generates a disk-patch in hydrostatic equilibrium
+# see Creasey, Theuns & Bower, 2013, for the equations:
+# disk parameters are: surface density sigma
+# scale height b
+# density: rho(z) = (sigma/2b) sech^2(z/b)
+# isothermal velocity dispersion = <v_z^2? = b pi G sigma
+# grad potential = 2 pi G sigma tanh(z/b)
+# potential = ln(cosh(z/b)) + const
+# Dynamical time = sqrt(b / (G sigma))
+# to obtain the 1/ch^2(z/b) profile from a uniform profile (a glass, say, or a uniform random variable), note that, when integrating in z
+# \int 0^z dz/ch^2(z) = tanh(z)-tanh(0) = \int_0^x dx = x (where the last integral refers to a uniform density distribution), so that z = atanh(x)
+# usage: python makeIC.py 1000
+
+# 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 = (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
+
+
+# parameters of potential
+surface_density = 10.
+scale_height = 400.
+gamma = 5./3.
+
+# 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
+utherm = math.pi * const_G * surface_density * scale_height / (gamma-1)
+v_disp = numpy.sqrt(2 * utherm)
+soundspeed = numpy.sqrt(utherm / (gamma * (gamma-1.)))
+t_dyn = numpy.sqrt(scale_height / (const_G * surface_density))
+t_cross = scale_height / soundspeed
+print 'dynamical time = ',t_dyn,' sound crossing time = ',t_cross,' sound speed= ',soundspeed,' 3D velocity dispersion = ',v_disp
+
+
+# Parameters
+periodic= 1 # 1 For periodic box
+boxSize = 400. # [kpc]
+Radius = 100. # maximum radius of particles [kpc]
+G = const_G
+
+# File
+fileName = "Disk-Patch.hdf5"
+
+#---------------------------------------------------
+mass = 1
+
+#--------------------------------------------------
+
+
+# read glass file and generate gas positions and tile it ntile times in each dimension
+inglass = '../../SodShock/glass_002.hdf5'
+infile = h5py.File(inglass, "r")
+one_glass_p = infile["/PartType0/Coordinates"][:,:]
+one_glass_h = infile["/PartType0/SmoothingLength"][:]
+ntile = 2
+
+# use fraction of these
+
+
+# scale in [-0.5,0.5]*BoxSize / ntile
+one_glass_p[:,:] -= 0.5
+one_glass_p *= boxSize / ntile
+one_glass_h *= boxSize / ntile
+ndens_glass = (one_glass_h.shape[0]) / (boxSize/ntile)**3
+h_glass = numpy.amin(one_glass_h) * (boxSize/ntile)
+
+#
+# glass_p = []
+# glass_h = []
+# for ix in range(0,ntile):
+# for iy in range(0,ntile):
+# for iz in range(0,ntile):
+# shift = one_glass_p.copy()
+# shift[:,0] += (ix-(ntile-1)/2.) * boxSize / ntile
+# shift[:,1] += (iy-(ntile-1)/2.) * boxSize / ntile
+# shift[:,2] += (iz-(ntile-1)/2.) * boxSize / ntile
+# glass_p.append(shift)
+# glass_h.append(one_glass_h.copy())
+
+# glass_p = numpy.concatenate(glass_p, axis=0)
+# glass_h = numpy.concatenate(glass_h, axis=0)
+
+# # generate density profile by rejecting particles based on density profile
+# prob = surface_density / (2.*scale_height) / numpy.cosh(abs(glass_p[:,2])/scale_height)**2
+# prob /= surface_density / (2.*scale_height) / numpy.cosh(0)**2
+# Nglass = glass_p.shape[0]
+# remain = 1. * numpy.random.rand(Nglass) < prob
+# for ic in range(3):
+# remain = remain & (glass_p[:,ic] > -boxSize/2) & (glass_p[:,ic] < boxSize/2)
+
+# #
+# numPart = 0
+# numGas = numpy.sum(remain)
+# pos = glass_p[remain,:]
+# h = glass_h[remain] * 1 / prob[remain]**(1./3.)
+# bad = h > boxSize/10.
+# h[bad] = boxSize/10.
+
+N = 1000
+Nran = 3*N
+r = numpy.zeros((N, 3))
+r[:,0] = (-1.+2*numpy.random.rand(N)) * boxSize / 2
+r[:,1] = (-1.+2*numpy.random.rand(N)) * boxSize / 2
+z = scale_height * numpy.arctanh(numpy.random.rand(Nran))
+gd = z < boxSize / 2
+r[:,2] = z[gd][0:N]
+random = numpy.random.rand(N) > 0.5
+r[random,2] *= -1
+pos = r.copy()
+r = 0
+numGas = numpy.shape(pos)[0]
+
+#
+column_density = surface_density * numpy.tanh(boxSize/2./scale_height)
+enclosed_mass = column_density * boxSize * boxSize
+pmass = enclosed_mass / numGas
+rho = surface_density / (2. * scale_height) / numpy.cosh(abs(pos[:,2])/scale_height)**2
+h = (pmass/rho)**(1./3.)
+bad = h > boxSize/10.
+h[bad] = boxSize/10.
+u = (1. + 0 * h) * utherm
+entropy = (gamma-1) * u / rho**(gamma-1)
+mass = 0.*h + pmass
+#mass = (1. + 0 * h) * surface_density / (2. * scale_height) / ndens_glass
+entropy_flag = 0
+vel = 0 + 0 * pos
+
+# move centre of disk to middle of box
+pos[:,:] += boxSize/2
+
+
+# create numPart dm particles
+numPart = 0
+
+# Create and write output file
+
+#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"] = [numGas, numPart, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numGas, 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"] = [entropy_flag]
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+
+# write gas particles
+grp0 = file.create_group("/PartType0")
+
+ds = grp0.create_dataset('Coordinates', (numGas, 3), 'f')
+ds[()] = pos
+
+ds = grp0.create_dataset('Velocities', (numGas, 3), 'f')
+ds[()] = vel
+
+ds = grp0.create_dataset('Masses', (numGas,), 'f')
+ds[()] = mass
+
+ds = grp0.create_dataset('SmoothingLength', (numGas,), 'f')
+ds[()] = h
+
+ds = grp0.create_dataset('InternalEnergy', (numGas,), 'f')
+u = numpy.full((numGas, ), utherm)
+if (entropy_flag == 1):
+ ds[()] = entropy
+else:
+ ds[()] = u
+
+ids = 1 + numpy.linspace(0, numGas, numGas, endpoint=False, dtype='L')
+ds = grp0.create_dataset('ParticleIDs', (numGas, ), 'L')
+ds[()] = ids
+
+# generate dark matter particles if needed
+if(numPart > 0):
+
+ # set seed for random number
+ numpy.random.seed(1234)
+
+#Particle group
+ 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))
+#
+ speed = vrot
+ v = numpy.zeros((numPart, 3))
+ omega = speed / radius
+ period = 2.*math.pi/omega
+ print 'period = minimum = ',min(period), ' maximum = ',max(period)
+
+ v[:,0] = -omega * r[:,1]
+ v[:,1] = omega * r[:,0]
+
+ ds = grp1.create_dataset('Coordinates', (numPart, 3), 'd')
+ ds[()] = r
+
+ ds = grp1.create_dataset('Velocities', (numPart, 3), 'f')
+ ds[()] = v
+ v = numpy.zeros(1)
+
+ m = numpy.full((numPart, ),10)
+ 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
+
+
+file.close()
+
+sys.exit()
diff --git a/examples/Disk-Patch/HydroStatic/makeIC.py b/examples/Disk-Patch/HydroStatic/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..cfc3ac638870535b7fb273f1e99bced2353e0f3f
--- /dev/null
+++ b/examples/Disk-Patch/HydroStatic/makeIC.py
@@ -0,0 +1,262 @@
+###############################################################################
+ # 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
+import matplotlib.pyplot as plt
+
+# Generates a disk-patch in hydrostatic equilibrium
+# see Creasey, Theuns & Bower, 2013, for the equations:
+# disk parameters are: surface density sigma
+# scale height b
+# density: rho(z) = (sigma/2b) sech^2(z/b)
+# isothermal velocity dispersion = <v_z^2? = b pi G sigma
+# grad potential = 2 pi G sigma tanh(z/b)
+# potential = ln(cosh(z/b)) + const
+# Dynamical time = sqrt(b / (G sigma))
+# to obtain the 1/ch^2(z/b) profile from a uniform profile (a glass, say, or a uniform random variable), note that, when integrating in z
+# \int 0^z dz/ch^2(z) = tanh(z)-tanh(0) = \int_0^x dx = x (where the last integral refers to a uniform density distribution), so that z = atanh(x)
+# usage: python makeIC.py 1000
+
+# 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 = (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
+
+
+# parameters of potential
+surface_density = 10.
+scale_height = 400.
+gamma = 5./3.
+
+# 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
+utherm = math.pi * const_G * surface_density * scale_height / (gamma-1)
+v_disp = numpy.sqrt(2 * utherm)
+soundspeed = numpy.sqrt(utherm / (gamma * (gamma-1.)))
+t_dyn = numpy.sqrt(scale_height / (const_G * surface_density))
+t_cross = scale_height / soundspeed
+print 'dynamical time = ',t_dyn,' sound crossing time = ',t_cross,' sound speed= ',soundspeed,' 3D velocity dispersion = ',v_disp
+
+
+# Parameters
+periodic= 1 # 1 For periodic box
+boxSize = 400. # [kpc]
+Radius = 100. # maximum radius of particles [kpc]
+G = const_G
+
+# File
+fileName = "Disk-Patch.hdf5"
+
+#---------------------------------------------------
+mass = 1
+
+#--------------------------------------------------
+
+
+# using glass ICs
+# read glass file and generate gas positions and tile it ntile times in each dimension
+ntile = 1
+inglass = '../../SodShock/glass_002.hdf5'
+infile = h5py.File(inglass, "r")
+one_glass_p = infile["/PartType0/Coordinates"][:,:]
+one_glass_h = infile["/PartType0/SmoothingLength"][:]
+
+# scale in [-0.5,0.5]*BoxSize / ntile
+one_glass_p[:,:] -= 0.5
+one_glass_p *= boxSize / ntile
+one_glass_h *= boxSize / ntile
+ndens_glass = (one_glass_h.shape[0]) / (boxSize/ntile)**3
+h_glass = numpy.amin(one_glass_h) * (boxSize/ntile)
+
+glass_p = []
+glass_h = []
+for ix in range(0,ntile):
+ for iy in range(0,ntile):
+ for iz in range(0,ntile):
+ shift = one_glass_p.copy()
+ shift[:,0] += (ix-(ntile-1)/2.) * boxSize / ntile
+ shift[:,1] += (iy-(ntile-1)/2.) * boxSize / ntile
+ shift[:,2] += (iz-(ntile-1)/2.) * boxSize / ntile
+ glass_p.append(shift)
+ glass_h.append(one_glass_h.copy())
+
+glass_p = numpy.concatenate(glass_p, axis=0)
+glass_h = numpy.concatenate(glass_h, axis=0)
+
+
+# use some of these
+numGas = 5000
+
+
+# use these as ICs
+pos = glass_p[0:numGas,:]
+h = glass_h[0:numGas]
+numGas = numpy.shape(pos)[0]
+print ' numGas = ', numGas
+
+
+# compute furthe properties of ICs
+column_density = surface_density * numpy.tanh(boxSize/2./scale_height)
+enclosed_mass = column_density * boxSize * boxSize
+pmass = enclosed_mass / numGas
+
+# desired density
+rho = surface_density / (2. * scale_height) / numpy.cosh(abs(pos[:,2])/scale_height)**2
+u = (1. + 0 * h) * utherm
+entropy = (gamma-1) * u / rho**(gamma-1)
+mass = 0.*h + pmass
+entropy_flag = 0
+vel = 0 + 0 * pos
+
+# move centre of disk to middle of box
+pos[:,:] += boxSize/2
+
+
+# create numPart dm particles
+numPart = 0
+
+# Create and write output file
+
+#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"] = [numGas, numPart, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numGas, 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"] = [entropy_flag]
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+
+# write gas particles
+grp0 = file.create_group("/PartType0")
+
+ds = grp0.create_dataset('Coordinates', (numGas, 3), 'f')
+ds[()] = pos
+
+ds = grp0.create_dataset('Velocities', (numGas, 3), 'f')
+ds[()] = vel
+
+ds = grp0.create_dataset('Masses', (numGas,), 'f')
+ds[()] = mass
+
+ds = grp0.create_dataset('SmoothingLength', (numGas,), 'f')
+ds[()] = h
+
+ds = grp0.create_dataset('InternalEnergy', (numGas,), 'f')
+u = numpy.full((numGas, ), utherm)
+if (entropy_flag == 1):
+ ds[()] = entropy
+else:
+ ds[()] = u
+
+ids = 1 + numpy.linspace(0, numGas, numGas, endpoint=False, dtype='L')
+ds = grp0.create_dataset('ParticleIDs', (numGas, ), 'L')
+ds[()] = ids
+
+# generate dark matter particles if needed
+if(numPart > 0):
+
+ # set seed for random number
+ numpy.random.seed(1234)
+
+#Particle group
+ 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))
+#
+ speed = vrot
+ v = numpy.zeros((numPart, 3))
+ omega = speed / radius
+ period = 2.*math.pi/omega
+ print 'period = minimum = ',min(period), ' maximum = ',max(period)
+
+ v[:,0] = -omega * r[:,1]
+ v[:,1] = omega * r[:,0]
+
+ ds = grp1.create_dataset('Coordinates', (numPart, 3), 'd')
+ ds[()] = r
+
+ ds = grp1.create_dataset('Velocities', (numPart, 3), 'f')
+ ds[()] = v
+ v = numpy.zeros(1)
+
+ m = numpy.full((numPart, ),10)
+ 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
+
+
+file.close()
+
+sys.exit()
diff --git a/examples/Disk-Patch/HydroStatic/test.pro b/examples/Disk-Patch/HydroStatic/test.pro
new file mode 100644
index 0000000000000000000000000000000000000000..675c30f6aabb1accabecb50f4314a1644c35d577
--- /dev/null
+++ b/examples/Disk-Patch/HydroStatic/test.pro
@@ -0,0 +1,193 @@
+;
+; 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 = 'Disk-Patch_'
+
+; set properties of potential
+uL = phys.pc ; unit of length
+uM = phys.msun ; unit of mass
+uV = 1d5 ; unit of velocity
+
+; properties of patch
+surface_density = 10.
+scale_height = 100.
+
+; derived units
+constG = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ;
+pcentre = [0.,0.,300.] * 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,'PartType0/ParticleIDs')
+nfollow = n_elements(id)
+
+; follow a subset
+nfollow = min(4000, nfollow) ; 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
+vzout = fltarr(nfollow, nsave) ; z
+rout = fltarr(nfollow, nsave) ; rho
+hout = fltarr(nfollow, nsave) ; h
+uout = fltarr(nfollow, nsave) ; thermal energy
+eout = fltarr(nfollow, nsave) ; energies
+ekin = fltarr(nfollow, nsave)
+epot = fltarr(nfollow, nsave) ; 2 pi G Sigma b ln(cosh(z/b)) + const
+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,'PartType0/Coordinates')
+ v = h5rd(inf,'PartType0/Velocities')
+ id = h5rd(inf,'PartType0/ParticleIDs')
+ rho = h5rd(inf,'PartType0/Density')
+ h = h5rd(inf,'PartType0/SmoothingLength')
+ utherm = h5rd(inf,'PartType0/InternalEnergy')
+ indx = sort(id)
+
+; if you want to sort particles by ID
+ id = id[indx]
+ rho = rho[indx]
+ utherm = utherm[indx]
+ h = h[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]
+ vzout[*,isave]= v[2,0:nfollow-1]
+ rout[*,isave] = rho[0:nfollow-1]
+ hout[*,isave] = h[0:nfollow-1]
+ uout[*,isave] = utherm[0:nfollow-1]
+ Lz = (p[0,*]-pcentre[0]) * v[1,*] - (p[1,*]-pcentre[1]) * v[0,*]
+ dz = reform(p[2,0:nfollow-1]-pcentre[2])
+; print,'time = ',time,p[0,0],v[0,0],id[0]
+ ek = 0.5 * dv
+ ep = fltarr(nfollow)
+ ep = 2 * !pi * constG * surface_density * scale_height * alog(cosh(abs(dz)/scale_height))
+ 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]
+ 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
+ nplot = nfollow
+
+ ; plot density profile
+ wset,0
+ xr = [0, 3*scale_height]
+ nbins = 100
+ zpos = findgen(nbins)/float(nbins-1) * max(xr)
+ dens = (surface_density/(2.d0*scale_height)) * 1./cosh(zpos/scale_height)^2
+ plot,[0],[0],xr=xr,/xs,yr=[0,max(dens)*1.4],/ys,/nodata,xtitle='|z|',ytitle='density'
+ oplot,zpos,dens,color=black,thick=3
+ oplot,abs(zout[*,1]),rout[*,1],psym=3
+ oplot,abs(zout[*,nsave-1]),rout[*,nsave-1],psym=3,color=red
+
+;; ; plot results on energy conservation for some particles
+;; nplot = min(100, 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 vertical oscillation
+ wset,2
+ xr = [min(tout), max(tout)]
+ yr = [-3,3]*scale_height
+ plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/nodata,xtitle='t',ytitle='z(t)'
+ color = floor(findgen(nplot)*255/float(nplot))
+ for i=0,nplot-1,50 do oplot,tout,zout[i,*],color=color(i)
+ screen_to_png,'orbit.png'
+
+
+; plot evolution of density for some particles close to disk
+ wset, 6
+ rr = (surface_density/(2.d0*scale_height)) * 1./cosh(abs(zout)/100.)^2 ; desired density
+ gd = where(abs(zout[*,0]) lt 50, ng)
+ plot,[0],[0],xr=[0, max(tout)], yr=10.^[-1,1], /xs, /ys, /nodata, /yl
+ nplot = min(40, ng)
+ color = floor(findgen(nplot)/float(nplot)*255)
+ for i=0, nplot-1 do oplot,tout[1:*],rout[gd[i],1:*]/rr[gd[i], 1:*], color=color[i],psym=-4
+
+
+;; ; 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
+
+
diff --git a/examples/Disk-Patch/configure b/examples/Disk-Patch/configure
new file mode 100644
index 0000000000000000000000000000000000000000..1cf691b1ebe3fcbf8667d501116e6c22a46b03e4
--- /dev/null
+++ b/examples/Disk-Patch/configure
@@ -0,0 +1 @@
+./configure --enable-optimization=no --enable-debug=yes --disable-vec
diff --git a/src/const.h b/src/const.h
index 34db85df64e023ca929b70e193747a3d472b98ff..c14419ac58f0130857643c76929cd56d5373a81a 100644
--- a/src/const.h
+++ b/src/const.h
@@ -68,6 +68,7 @@
//#define EXTERNAL_POTENTIAL_POINTMASS
//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
#define EXTERNAL_POTENTIAL_DISK_PATCH
+#define ISOTHERMAL_GLASS
/* Are we debugging ? */
//#define SWIFT_DEBUG_CHECKS
diff --git a/src/kick.h b/src/kick.h
index b57e13d4ebf27d3a366d571e7fd4cd819653f726..7c6dbd16fbf613b63ca3225a03ebb73f2a486b36 100644
--- a/src/kick.h
+++ b/src/kick.h
@@ -46,11 +46,31 @@ __attribute__((always_inline)) INLINE static void kick_gpart(
gp->ti_begin = gp->ti_end;
gp->ti_end = gp->ti_begin + new_dti;
+
+#ifdef ISOTHERMAL_GLASS
+ /* TT: add viscous drag */
+ float a_tot[3] = {gp->a_grav[0], gp->a_grav[1], gp->a_grav[2]};
+ const float surface_density= 10.;
+ const float scale_height = 100;
+ const float t_dyn = sqrt(scale_height / (const_G * surface_density));
+ const double time = ti_start * timeBase;
+ a_tot[0] -= gp->v_full[0] / (0.05 * t_dyn);
+ a_tot[1] -= gp->v_full[1] / (0.05 * t_dyn);
+ a_tot[2] -= gp->v_full[2] / (0.05 * t_dyn);
+
+ /* Kick particles in momentum space */
+ gp->v_full[0] += a_tot[0] * dt;
+ gp->v_full[1] += a_tot[1] * dt;
+ gp->v_full[2] += a_tot[2] * dt;
+
+#else
+
/* Kick particles in momentum space */
gp->v_full[0] += gp->a_grav[0] * dt;
gp->v_full[1] += gp->a_grav[1] * dt;
gp->v_full[2] += gp->a_grav[2] * dt;
+#endif
/* Extra kick work */
gravity_kick_extra(gp, dt, half_dt);
}
@@ -89,6 +109,17 @@ __attribute__((always_inline)) INLINE static void kick_part(
a_tot[2] += p->gpart->a_grav[2];
}
+#ifdef ISOTHERMAL_GLASS
+ /* TT: add viscous drag */
+ const float surface_density= 10.;
+ const float scale_height = 100;
+ const float t_dyn = sqrt(scale_height / (const_G * surface_density));
+ const double time = ti_start * timeBase;
+ a_tot[0] -= p->gpart->v_full[0] / (0.05 * t_dyn);
+ a_tot[1] -= p->gpart->v_full[1] / (0.05 * t_dyn);
+ a_tot[2] -= p->gpart->v_full[2] / (0.05 * t_dyn);
+#endif
+
/* Kick particles in momentum space */
xp->v_full[0] += a_tot[0] * dt;
xp->v_full[1] += a_tot[1] * dt;
diff --git a/src/potentials.h b/src/potentials.h
index 83bdaf6fc817344f97612a8e8ddfc42c2d7ffcfa..cc49545e95b9a8859277028849bd86b73ce03391 100644
--- a/src/potentials.h
+++ b/src/potentials.h
@@ -26,6 +26,7 @@
/* Some standard headers. */
#include <math.h>
+#include <float.h>
/* Local includes. */
#include "const.h"
@@ -35,7 +36,8 @@
#include "physical_constants.h"
#include "units.h"
#include "math.h"
-
+#define FINLINE __attribute__((always_inline)) INLINE
+/* #define FINLINE */
/* External Potential Properties */
struct external_potential {
@@ -72,27 +74,55 @@ struct external_potential {
* @param phys_cont The physical constants in internal units.
* @param gp Pointer to the g-particle data.
*/
-__attribute__((always_inline))
- INLINE static float external_gravity_disk_patch_timestep(
+FINLINE static float external_gravity_disk_patch_timestep(
const struct external_potential* potential,
const struct phys_const* const phys_const,
const struct gpart* const g) {
+
+ /* initilize time step to disk dynamical time */
+ const float dt_dyn = sqrt(potential->disk_patch_potential.scale_height / (phys_const->const_newton_G * potential->disk_patch_potential.surface_density));
+ float dt = dt_dyn;
+
+
+ /* absolute value of height above disk */
const float dz = fabs(g->x[2] - potential->disk_patch_potential.z_disk);
+
+ /* vertical cceleration */
const float z_accel = 2 * M_PI * phys_const->const_newton_G * potential->disk_patch_potential.surface_density
- * tanh(dz / potential->disk_patch_potential.scale_height);
- const float dz_accel_over_dt = 2 * M_PI * phys_const->const_newton_G * potential->disk_patch_potential.surface_density /
- potential->disk_patch_potential.scale_height / cosh(dz / potential->disk_patch_potential.scale_height) / cosh(dz / potential->disk_patch_potential.scale_height) * fabs(g->v_full[2]);
+ * tanh(dz / potential->disk_patch_potential.scale_height);
+
+ /* demand that dt * velocity < fraction of scale height of disk */
+ float dt1 = FLT_MAX;
+ if(fabs(g->v_full[2]) > 0)
+ {
+ dt1 = potential->disk_patch_potential.scale_height / fabs(g->v_full[2]);
+ if(dt1 < dt) dt = dt1;
+ }
- /* demand that time step * derivative of acceleration < fraction of acceleration */
- const float dt1 = potential->disk_patch_potential.timestep_mult * z_accel / dz_accel_over_dt;
-
/* demand that dt^2 * acceleration < fraction of scale height of disk */
- const float dt2 = potential->disk_patch_potential.timestep_mult * potential->disk_patch_potential.scale_height / fabs(z_accel);
-
- float dt = dt1;
- if(dt2 < dt) dt = dt2;
- return dt2;
+ float dt2 = FLT_MAX;
+ if(fabs(z_accel) > 0)
+ {
+ dt2 = potential->disk_patch_potential.scale_height / fabs(z_accel);
+ if(dt2 < dt * dt) dt = sqrt(dt2);
+ }
+
+ /* demand that dt^3 jerk < fraction of scale height of disk */
+ float dt3 = FLT_MAX;
+ if(abs(g->v_full[2]) > 0)
+ {
+ const float dz_accel_over_dt = 2 * M_PI * phys_const->const_newton_G * potential->disk_patch_potential.surface_density /
+ potential->disk_patch_potential.scale_height / cosh(dz / potential->disk_patch_potential.scale_height) / cosh(dz / potential->disk_patch_potential.scale_height) * fabs(g->v_full[2]);
+
+ dt3 = potential->disk_patch_potential.scale_height / fabs(dz_accel_over_dt);
+ if(dt3 < dt * dt * dt) dt = pow(dt3, 1./3.);
+ }
+
+ // if(potential->disk_patch_potential.timestep_mult * dt < 1e-3 * dt_dyn)
+ // message(" i= %lld t_dyn= %e dt_vel= %e dt_accel= %e dt_jerk= %e", g->id_or_neg_offset, dt_dyn, dt1, dt2, dt3);
+
+ return potential->disk_patch_potential.timestep_mult * dt;
}
/**
* @brief Computes the gravitational acceleration from a hydrostatic disk (Creasy, Theuns & Bower, 2013)
@@ -100,13 +130,17 @@ __attribute__((always_inline))
* @param phys_cont The physical constants in internal units.
* @param gp Pointer to the g-particle data.
*/
-__attribute__((always_inline)) INLINE static void external_gravity_disk_patch_potential(const struct external_potential* potential, const struct phys_const* const phys_const, struct gpart* g) {
+FINLINE static void external_gravity_disk_patch_potential(const struct external_potential* potential, const struct phys_const* const phys_const, struct gpart* g) {
const float dz = g->x[2] - potential->disk_patch_potential.z_disk;
- g->a_grav[0] = 0;
- g->a_grav[1] = 0;
- const float z_accel = 2 * M_PI * phys_const->const_newton_G * potential->disk_patch_potential.surface_density
- * tanh(fabs(dz) / potential->disk_patch_potential.scale_height);
+ g->a_grav[0] += 0;
+ g->a_grav[1] += 0;
+
+ /* const float z_accel = 2 * M_PI * phys_const->const_newton_G * potential->disk_patch_potential.surface_density */
+ /* * tanh(fabs(dz) / potential->disk_patch_potential.scale_height); */
+ /* impose *no* graitational acceleration */
+ const float z_accel = 0.;
+
if (dz > 0)
g->a_grav[2] -= z_accel;
if (dz < 0)
@@ -127,6 +161,7 @@ __attribute__((always_inline)) INLINE static float
external_gravity_isothermalpotential_timestep(
const struct external_potential* potential,
const struct phys_const* const phys_const, const struct gpart* const g) {
+ const struct gpart* const g) {
const float dx = g->x[0] - potential->isothermal_potential.x;
const float dy = g->x[1] - potential->isothermal_potential.y;
@@ -192,9 +227,8 @@ external_gravity_isothermalpotential(const struct external_potential* potential,
* @param phys_const The physical constants in internal units.
* @param g Pointer to the g-particle data.
*/
-__attribute__((always_inline)) INLINE static float
-external_gravity_pointmass_timestep(const struct external_potential* potential,
- const struct phys_const* const phys_const,
+FINLINE static float external_gravity_pointmass_timestep(const struct external_potential* potential,
+ const struct phys_const* const phys_const,
const struct gpart* const g) {
const float G_newton = phys_const->const_newton_G;
@@ -228,7 +262,7 @@ external_gravity_pointmass_timestep(const struct external_potential* potential,
* @param phys_const The physical constants in internal units.
* @param g Pointer to the g-particle data.
*/
-__attribute__((always_inline)) INLINE static void external_gravity_pointmass(
+FINLINE static void external_gravity_pointmass(
const struct external_potential* potential,
const struct phys_const* const phys_const, struct gpart* g) {