added initial conditions for disk-patch

examples/Disk-Patch/GravityOnly/disk-patch.yml

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

examples/Disk-Patch/GravityOnly/legend.pro

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

examples/Disk-Patch/GravityOnly/makeIC.py

+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ #                    Tom Theuns (tom.theuns@durham.ac.uk)
+ # 
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ # 
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ # GNU General Public License for more details.
+ # 
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ # 
+ ##############################################################################
+
+import h5py
+import sys
+import numpy
+import math
+import random
+
+# Generates N particles in a 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()

examples/Disk-Patch/GravityOnly/readme.txt

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

examples/Disk-Patch/GravityOnly/run.sh

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

examples/Disk-Patch/GravityOnly/test.pro

+;
+;  test energy / angular momentum conservation of test problem
+;
+
+iplot = 1 ; if iplot = 1, make plot of E/Lz conservation, else, simply compare final and initial energy
+
+; set physical constants
+@physunits
+
+indir    = './'
+basefile = '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
+
+

examples/Disk-Patch/HydroStatic/disk-patch.yml

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

examples/Disk-Patch/HydroStatic/legend.pro

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

examples/Disk-Patch/HydroStatic/makeIC-stretch.py

+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ #                    Tom Theuns (tom.theuns@durham.ac.uk)
+ # 
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ # 
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ # GNU General Public License for more details.
+ # 
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ # 
+ ##############################################################################
+
+import h5py
+import sys
+import numpy
+import math
+import random
+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()

examples/Disk-Patch/HydroStatic/makeIC.py

+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ #                    Tom Theuns (tom.theuns@durham.ac.uk)
+ # 
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ # 
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ # GNU General Public License for more details.
+ # 
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ # 
+ ##############################################################################
+
+import h5py
+import sys
+import numpy
+import math
+import random
+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()

examples/Disk-Patch/HydroStatic/test.pro

+;
+;  test energy / angular momentum conservation of test problem
+;
+
+iplot = 1 ; if iplot = 1, make plot of E/Lz conservation, else, simply compare final and initial energy
+
+; set physical constants
+@physunits
+
+indir    = './'
+basefile = '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
+
+

examples/Disk-Patch/configure

+./configure --enable-optimization=no --enable-debug=yes --disable-vec

src/const.h

@@ -66,6 +66,7 @@
 /* #define EXTERNAL_POTENTIAL_POINTMASS */
 //#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
 #define EXTERNAL_POTENTIAL_DISK_PATCH
+#define ISOTHERMAL_GLASS
 
 /* Are we debugging ? */
 //#define SWIFT_DEBUG_CHECKS

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);
 }
@@ -86,9 +106,20 @@ __attribute__((always_inline)) INLINE static void kick_part(
   if (p->gpart != NULL) {
     a_tot[0] += p->gpart->a_grav[0];
     a_tot[1] += p->gpart->a_grav[1];
-    a_tot[1] += p->gpart->a_grav[2];
+    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;

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 {
 
@@ -71,27 +73,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)
@@ -99,13 +129,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)
@@ -122,8 +156,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_disk_patch_po
  * @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_isothermalpotential_timestep(
+FINLINE static float external_gravity_isothermalpotential_timestep(
 		  const struct external_potential* potential,
         const struct phys_const* const phys_const,
         const struct gpart* const g) {
@@ -154,7 +187,7 @@ __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_isothermalpotential(const struct external_potential* potential, const struct phys_const* const phys_const, struct gpart* g) {
+FINLINE static void external_gravity_isothermalpotential(const struct external_potential* potential, const struct phys_const* const phys_const, struct gpart* g) {
 
   const float dx    = g->x[0] - potential->isothermal_potential.x;
   const float dy    = g->x[1] - potential->isothermal_potential.y;
@@ -181,9 +214,8 @@ __attribute__((always_inline)) INLINE static void external_gravity_isothermalpot
  * @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;
@@ -216,7 +248,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) {