diff --git a/examples/IsothermalPotential/GravityOnly/externalGravity.yml b/examples/IsothermalPotential/GravityOnly/externalGravity.yml
new file mode 100644
index 0000000000000000000000000000000000000000..53c8a0cf20a0fd9666b5f9fb43baeff6045ce978
--- /dev/null
+++ b/examples/IsothermalPotential/GravityOnly/externalGravity.yml
@@ -0,0 +1,73 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.9885e33 # Grams
+ UnitLength_in_cgs: 3.0856776e21 # Centimeters
+ UnitVelocity_in_cgs: 1e5 # Centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 8. # The end time of the simulation (in internal units).
+ dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-1 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-2 # Time between statistics output
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: Isothermal # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 0.02 # Time difference between consecutive outputs (in internal units)
+
+# Parameters for the task scheduling
+Scheduler:
+ nr_threads: 8 # The number of threads per MPI rank to use.
+ nr_queues: 0 # The number of task queues to use. Use 0 to let the system decide.
+ cell_max_size: 8000000 # Maximal number of interactions per task (this is the default value).
+ cell_sub_size: 8000000 # Maximal number of interactions per sub-task (this is the default value).
+ cell_split_size: 400 # Maximal number of particles per cell (this is the default value).
+
+
+# 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: Isothermal.hdf5 # The file to read
+ h_scaling: 1. # A scaling factor to apply to all smoothing lengths in the ICs.
+ shift_x: 100. # A shift to apply to all particles read from the ICs (in internal units).
+ shift_y: 100.
+ shift_z: 100.
+
+# Parameters govering domain decomposition
+DomainDecomposition:
+ initial_type: m # The initial strategy ("g", "m", "w", or "v"). See documentation for details.
+ initial_grid_x: 10 # Grid size if the 'g' strategy is chosen.
+ initial_grid_y: 10
+ initial_grid_z: 10
+ repartition_type: b # The re-decomposition strategy ("n", "b", "v", "e" or "x"). See documentation for details.
+
+
+# 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
+
+IsothermalPotential:
+ position_x: 100. # location of external point mass in internal units
+ position_y: 100.
+ position_z: 100.
+ vrot: 200. # rotation speed of isothermal potential in internal units
+
+
diff --git a/examples/IsothermalPotential/GravityOnly/legend.pro b/examples/IsothermalPotential/GravityOnly/legend.pro
new file mode 100755
index 0000000000000000000000000000000000000000..e510bdfcd7dcabed796cb090c6fa7d54536608b6
--- /dev/null
+++ b/examples/IsothermalPotential/GravityOnly/legend.pro
@@ -0,0 +1,459 @@
+;+
+; NAME:
+; LEGEND
+; PURPOSE:
+; Create an annotation legend for a plot.
+; EXPLANATION:
+; This procedure makes a legend for a plot. The legend can contain
+; a mixture of symbols, linestyles, Hershey characters (vectorfont),
+; and filled polygons (usersym). A test procedure, legendtest.pro,
+; shows legend's capabilities. Placement of the legend is controlled
+; with keywords like /right, /top, and /center or by using a position
+; keyword for exact placement (position=[x,y]) or via mouse (/position).
+; CALLING SEQUENCE:
+; LEGEND [,items][,keyword options]
+; EXAMPLES:
+; The call:
+; legend,['Plus sign','Asterisk','Period'],psym=[1,2,3]
+; produces:
+; -----------------
+; | |
+; | + Plus sign |
+; | * Asterisk |
+; | . Period |
+; | |
+; -----------------
+; Each symbol is drawn with a plots command, so they look OK.
+; Other examples are given in optional output keywords.
+;
+; lines = indgen(6) ; for line styles
+; items = 'linestyle '+strtrim(lines,2) ; annotations
+; legend,items,linestyle=lines ; vertical legend---upper left
+; items = ['Plus sign','Asterisk','Period']
+; sym = [1,2,3]
+; legend,items,psym=sym ; ditto except using symbols
+; legend,items,psym=sym,/horizontal ; horizontal format
+; legend,items,psym=sym,box=0 ; sans border
+; legend,items,psym=sym,delimiter='=' ; embed '=' betw psym & text
+; legend,items,psym=sym,margin=2 ; 2-character margin
+; legend,items,psym=sym,position=[x,y] ; upper left in data coords
+; legend,items,psym=sym,pos=[x,y],/norm ; upper left in normal coords
+; legend,items,psym=sym,pos=[x,y],/device ; upper left in device coords
+; legend,items,psym=sym,/position ; interactive position
+; legend,items,psym=sym,/right ; at upper right
+; legend,items,psym=sym,/bottom ; at lower left
+; legend,items,psym=sym,/center ; approximately near center
+; legend,items,psym=sym,number=2 ; plot two symbols, not one
+; legend,items,/fill,psym=[8,8,8],colors=[10,20,30]; 3 filled squares
+; INPUTS:
+; items = text for the items in the legend, a string array.
+; For example, items = ['diamond','asterisk','square'].
+; You can omit items if you don't want any text labels.
+; OPTIONAL INPUT KEYWORDS:
+;
+; linestyle = array of linestyle numbers If linestyle[i] < 0, then omit
+; ith symbol or line to allow a multi-line entry. If
+; linestyle = -99 then text will be left-justified.
+; psym = array of plot symbol numbers. If psym[i] is negative, then a
+; line connects pts for ith item. If psym[i] = 8, then the
+; procedure usersym is called with vertices define in the
+; keyword usersym. If psym[i] = 88, then use the previously
+; defined user symbol
+; vectorfont = vector-drawn characters for the sym/line column, e.g.,
+; ['!9B!3','!9C!3','!9D!3'] produces an open square, a checkmark,
+; and a partial derivative, which might have accompanying items
+; ['BOX','CHECK','PARTIAL DERIVATIVE'].
+; There is no check that !p.font is set properly, e.g., -1 for
+; X and 0 for PostScript. This can produce an error, e.g., use
+; !20 with PostScript and !p.font=0, but allows use of Hershey
+; *AND* PostScript fonts together.
+; N. B.: Choose any of linestyle, psym, and/or vectorfont. If none is
+; present, only the text is output. If more than one
+; is present, all need the same number of elements, and normal
+; plot behaviour occurs.
+; By default, if psym is positive, you get one point so there is
+; no connecting line. If vectorfont[i] = '',
+; then plots is called to make a symbol or a line, but if
+; vectorfont[i] is a non-null string, then xyouts is called.
+; /help = flag to print header
+; /horizontal = flag to make the legend horizontal
+; /vertical = flag to make the legend vertical (D=vertical)
+; box = flag to include/omit box around the legend (D=include)
+; clear = flag to clear the box area before drawing the legend
+; delimiter = embedded character(s) between symbol and text (D=none)
+; colors = array of colors for plot symbols/lines (D=!P.color)
+; textcolors = array of colors for text (D=!P.color)
+; margin = margin around text measured in characters and lines
+; spacing = line spacing (D=bit more than character height)
+; pspacing = psym spacing (D=3 characters) (when number of symbols is
+; greater than 1)
+; charsize = just like !p.charsize for plot labels
+; charthick = just like !p.charthick for plot labels
+; thick = array of line thickness numbers (D = !P.thick), if used, then
+; linestyle must also be specified
+; position = data coordinates of the /top (D) /left (D) of the legend
+; normal = use normal coordinates for position, not data
+; device = use device coordinates for position, not data
+; number = number of plot symbols to plot or length of line (D=1)
+; usersym = 2-D array of vertices, cf. usersym in IDL manual.
+; (/USERSYM =square, default is to use existing USERSYM definition)
+; /fill = flag to fill the usersym
+; /left_legend = flag to place legend snug against left side of plot
+; window (D)
+; /right_legend = flag to place legend snug against right side of plot
+; window. If /right,pos=[x,y], then x is position of RHS and
+; text runs right-to-left.
+; /top_legend = flag to place legend snug against top of plot window (D)
+; /bottom = flag to place legend snug against bottom of plot window
+; /top,pos=[x,y] and /bottom,pos=[x,y] produce same positions.
+;
+; If LINESTYLE, PSYM, VECTORFONT, THICK, COLORS, or TEXTCOLORS are
+; supplied as scalars, then the scalar value is set for every line or
+; symbol in the legend.
+; Outputs:
+; legend to current plot device
+; OPTIONAL OUTPUT KEYWORDS:
+; corners = 4-element array, like !p.position, of the normalized
+; coords for the box (even if box=0): [llx,lly,urx,ury].
+; Useful for multi-column or multi-line legends, for example,
+; to make a 2-column legend, you might do the following:
+; c1_items = ['diamond','asterisk','square']
+; c1_psym = [4,2,6]
+; c2_items = ['solid','dashed','dotted']
+; c2_line = [0,2,1]
+; legend,c1_items,psym=c1_psym,corners=c1,box=0
+; legend,c2_items,line=c2_line,corners=c2,box=0,pos=[c1[2],c1[3]]
+; c = [c1[0]<c2[0],c1[1]<c2[1],c1[2]>c2[2],c1[3]>c2[3]]
+; plots,[c[0],c[0],c[2],c[2],c[0]],[c[1],c[3],c[3],c[1],c[1]],/norm
+; Useful also to place the legend. Here's an automatic way to place
+; the legend in the lower right corner. The difficulty is that the
+; legend's width is unknown until it is plotted. In this example,
+; the legend is plotted twice: the first time in the upper left, the
+; second time in the lower right.
+; legend,['1','22','333','4444'],linestyle=indgen(4),corners=corners
+; ; BOGUS LEGEND---FIRST TIME TO REPORT CORNERS
+; xydims = [corners[2]-corners[0],corners[3]-corners[1]]
+; ; SAVE WIDTH AND HEIGHT
+; chdim=[!d.x_ch_size/float(!d.x_size),!d.y_ch_size/float(!d.y_size)]
+; ; DIMENSIONS OF ONE CHARACTER IN NORMALIZED COORDS
+; pos = [!x.window[1]-chdim[0]-xydims[0] $
+; ,!y.window[0]+chdim[1]+xydims[1]]
+; ; CALCULATE POSITION FOR LOWER RIGHT
+; plot,findgen(10) ; SIMPLE PLOT; YOU DO WHATEVER YOU WANT HERE.
+; legend,['1','22','333','4444'],linestyle=indgen(4),pos=pos
+; ; REDO THE LEGEND IN LOWER RIGHT CORNER
+; You can modify the pos calculation to place the legend where you
+; want. For example to place it in the upper right:
+; pos = [!x.window[1]-chdim[0]-xydims[0],!y.window[1]-xydims[1]]
+; Common blocks:
+; none
+; Procedure:
+; If keyword help is set, call doc_library to print header.
+; See notes in the code. Much of the code deals with placement of the
+; legend. The main problem with placement is not being
+; able to sense the length of a string before it is output. Some crude
+; approximations are used for centering.
+; Restrictions:
+; Here are some things that aren't implemented.
+; - An orientation keyword would allow lines at angles in the legend.
+; - An array of usersyms would be nice---simple change.
+; - An order option to interchange symbols and text might be nice.
+; - Somebody might like double boxes, e.g., with box = 2.
+; - Another feature might be a continuous bar with ticks and text.
+; - There are no guards to avoid writing outside the plot area.
+; - There is no provision for multi-line text, e.g., '1st line!c2nd line'
+; Sensing !c would be easy, but !c isn't implemented for PostScript.
+; A better way might be to simply output the 2nd line as another item
+; but without any accompanying symbol or linestyle. A flag to omit
+; the symbol and linestyle is linestyle[i] = -1.
+; - There is no ability to make a title line containing any of titles
+; for the legend, for the symbols, or for the text.
+; Side Effects:
+; Modification history:
+; write, 24-25 Aug 92, F K Knight (knight@ll.mit.edu)
+; allow omission of items or omission of both psym and linestyle, add
+; corners keyword to facilitate multi-column legends, improve place-
+; ment of symbols and text, add guards for unequal size, 26 Aug 92, FKK
+; add linestyle(i)=-1 to suppress a single symbol/line, 27 Aug 92, FKK
+; add keyword vectorfont to allow characters in the sym/line column,
+; 28 Aug 92, FKK
+; add /top, /bottom, /left, /right keywords for automatic placement at
+; the four corners of the plot window. The /right keyword forces
+; right-to-left printing of menu. 18 Jun 93, FKK
+; change default position to data coords and add normal, data, and
+; device keywords, 17 Jan 94, FKK
+; add /center keyword for positioning, but it is not precise because
+; text string lengths cannot be known in advance, 17 Jan 94, FKK
+; add interactive positioning with /position keyword, 17 Jan 94, FKK
+; allow a legend with just text, no plotting symbols. This helps in
+; simply describing a plot or writing assumptions done, 4 Feb 94, FKK
+; added thick, symsize, and clear keyword Feb 96, W. Landsman HSTX
+; David Seed, HR Wallingford, d.seed@hrwallingford.co.uk
+; allow scalar specification of keywords, Mar 96, W. Landsman HSTX
+; added charthick keyword, June 96, W. Landsman HSTX
+; Made keyword names left,right,top,bottom,center longer,
+; Aug 16, 2000, Kim Tolbert
+; Added ability to have regular text lines in addition to plot legend
+; lines in legend. If linestyle is -99 that item is left-justified.
+; Previously, only option for no sym/line was linestyle=-1, but then text
+; was lined up after sym/line column. 10 Oct 2000, Kim Tolbert
+; Make default value of thick = !P.thick W. Landsman Jan. 2001
+; Don't overwrite existing USERSYM definition W. Landsman Mar. 2002
+;-
+pro legend, items, BOTTOM_LEGEND=bottom, BOX = box, CENTER_LEGEND=center, $
+ CHARTHICK=charthick, CHARSIZE = charsize, CLEAR = clear, COLORS = colorsi, $
+ CORNERS = corners, DATA=data, DELIMITER=delimiter, DEVICE=device, $
+ FILL=fill, HELP = help, HORIZONTAL=horizontal,LEFT_LEGEND=left, $
+ LINESTYLE=linestylei, MARGIN=margin, NORMAL=normal, NUMBER=number, $
+ POSITION=position,PSPACING=pspacing, PSYM=psymi, RIGHT_LEGEND=right, $
+ SPACING=spacing, SYMSIZE=symsize, TEXTCOLORS=textcolorsi, THICK=thicki, $
+ TOP_LEGEND=top, USERSYM=usersym, VECTORFONT=vectorfonti, VERTICAL=vertical
+;
+; =====>> HELP
+;
+on_error,2
+if keyword_set(help) then begin & doc_library,'legend' & return & endif
+;
+; =====>> SET DEFAULTS FOR SYMBOLS, LINESTYLES, AND ITEMS.
+;
+ ni = n_elements(items)
+ np = n_elements(psymi)
+ nl = n_elements(linestylei)
+ nth = n_elements(thicki)
+ nv = n_elements(vectorfonti)
+ nlpv = max([np,nl,nv])
+ n = max([ni,np,nl,nv]) ; NUMBER OF ENTRIES
+strn = strtrim(n,2) ; FOR ERROR MESSAGES
+if n eq 0 then message,'No inputs! For help, type legend,/help.'
+if ni eq 0 then begin
+ items = replicate('',n) ; DEFAULT BLANK ARRAY
+endif else begin
+ if size(items,/TNAME) NE 'STRING' then message, $
+ 'First parameter must be a string array. For help, type legend,/help.'
+ if ni ne n then message,'Must have number of items equal to '+strn
+endelse
+symline = (np ne 0) or (nl ne 0) ; FLAG TO PLOT SYM/LINE
+ if (np ne 0) and (np ne n) and (np NE 1) then message, $
+ 'Must have 0, 1 or '+strn+' elements in PSYM array.'
+ if (nl ne 0) and (nl ne n) and (nl NE 1) then message, $
+ 'Must have 0, 1 or '+strn+' elements in LINESTYLE array.'
+ if (nth ne 0) and (nth ne n) and (nth NE 1) then message, $
+ 'Must have 0, 1 or '+strn+' elements in THICK array.'
+
+ case nl of
+ 0: linestyle = intarr(n) ;Default = solid
+ 1: linestyle = intarr(n) + linestylei
+ else: linestyle = linestylei
+ endcase
+
+ case nth of
+ 0: thick = replicate(!p.thick,n) ;Default = !P.THICK
+ 1: thick = intarr(n) + thicki
+ else: thick = thicki
+ endcase
+
+ case np of ;Get symbols
+ 0: psym = intarr(n) ;Default = solid
+ 1: psym = intarr(n) + psymi
+ else: psym = psymi
+ endcase
+
+ case nv of
+ 0: vectorfont = replicate('',n)
+ 1: vectorfont = replicate(vectorfonti,n)
+ else: vectorfont = vectorfonti
+ endcase
+;
+; =====>> CHOOSE VERTICAL OR HORIZONTAL ORIENTATION.
+;
+if n_elements(horizontal) eq 0 then begin ; D=VERTICAL
+ if n_elements(vertical) eq 0 then vertical = 1
+endif else begin
+ if n_elements(vertical) eq 0 then vertical = not horizontal
+endelse
+;
+; =====>> SET DEFAULTS FOR OTHER OPTIONS.
+;
+if n_elements(box) eq 0 then box = 1
+if n_elements(clear) eq 0 then clear = 0
+
+if n_elements(margin) eq 0 then margin = 0.5
+if n_elements(delimiter) eq 0 then delimiter = ''
+if n_elements(charsize) eq 0 then charsize = !p.charsize
+if n_elements(charthick) eq 0 then charthick = !p.charthick
+if charsize eq 0 then charsize = 1
+if (n_elements (symsize) eq 0) then symsize= charsize + intarr(n)
+if n_elements(number) eq 0 then number = 1
+ case N_elements(colorsi) of
+ 0: colors = replicate(!P.color,n) ;Default is !P.COLOR
+ 1: colors = replicate(colorsi,n)
+ else: colors = colorsi
+ endcase
+
+ case N_elements(textcolorsi) of
+ 0: textcolors = replicate(!P.color,n) ;Default is !P.COLOR
+ 1: textcolors = replicate(textcolorsi,n)
+ else: textcolors = textcolorsi
+ endcase
+ fill = keyword_set(fill)
+if n_elements(usersym) eq 1 then usersym = 2*[[0,0],[0,1],[1,1],[1,0],[0,0]]-1
+;
+; =====>> INITIALIZE SPACING
+;
+if n_elements(spacing) eq 0 then spacing = 1.2
+if n_elements(pspacing) eq 0 then pspacing = 3
+xspacing = !d.x_ch_size/float(!d.x_size) * (spacing > charsize)
+yspacing = !d.y_ch_size/float(!d.y_size) * (spacing > charsize)
+ltor = 1 ; flag for left-to-right
+if n_elements(left) eq 1 then ltor = left eq 1
+if n_elements(right) eq 1 then ltor = right ne 1
+ttob = 1 ; flag for top-to-bottom
+if n_elements(top) eq 1 then ttob = top eq 1
+if n_elements(bottom) eq 1 then ttob = bottom ne 1
+xalign = ltor ne 1 ; x alignment: 1 or 0
+yalign = -0.5*ttob + 1 ; y alignment: 0.5 or 1
+xsign = 2*ltor - 1 ; xspacing direction: 1 or -1
+ysign = 2*ttob - 1 ; yspacing direction: 1 or -1
+if not ttob then yspacing = -yspacing
+if not ltor then xspacing = -xspacing
+;
+; =====>> INITIALIZE POSITIONS: FIRST CALCULATE X OFFSET FOR TEXT
+;
+xt = 0
+if nlpv gt 0 then begin ; SKIP IF TEXT ITEMS ONLY.
+if vertical then begin ; CALC OFFSET FOR TEXT START
+ for i = 0,n-1 do begin
+ if (psym[i] eq 0) and (vectorfont[i] eq '') then num = (number + 1) > 3 else num = number
+ if psym[i] lt 0 then num = number > 2 ; TO SHOW CONNECTING LINE
+ if psym[i] eq 0 then expand = 1 else expand = 2
+ thisxt = (expand*pspacing*(num-1)*xspacing)
+ if ltor then xt = thisxt > xt else xt = thisxt < xt
+ endfor
+endif ; NOW xt IS AN X OFFSET TO ALIGN ALL TEXT ENTRIES.
+endif
+;
+; =====>> INITIALIZE POSITIONS: SECOND LOCATE BORDER
+;
+if !x.window[0] eq !x.window[1] then begin
+ plot,/nodata,xstyle=4,ystyle=4,[0],/noerase
+endif
+; next line takes care of weirdness with small windows
+pos = [min(!x.window),min(!y.window),max(!x.window),max(!y.window)]
+case n_elements(position) of
+ 0: begin
+ if ltor then px = pos[0] else px = pos[2]
+ if ttob then py = pos[3] else py = pos[1]
+ if keyword_set(center) then begin
+ if not keyword_set(right) and not keyword_set(left) then $
+ px = (pos[0] + pos[2])/2. - xt
+ if not keyword_set(top) and not keyword_set(bottom) then $
+ py = (pos[1] + pos[3])/2. + n*yspacing
+ endif
+ position = [px,py] + [xspacing,-yspacing]
+ end
+ 1: begin ; interactive
+ message,/inform,'Place mouse at upper left corner and click any mouse button.'
+ cursor,x,y,/normal
+ position = [x,y]
+ end
+ 2: begin ; convert upper left corner to normal coordinates
+ if keyword_set(data) then $
+ position = convert_coord(position,/to_norm) $
+ else if keyword_set(device) then $
+ position = convert_coord(position,/to_norm,/device) $
+ else if not keyword_set(normal) then $
+ position = convert_coord(position,/to_norm)
+ end
+ else: message,'Position keyword can have 0, 1, or 2 elements only. Try legend,/help.'
+endcase
+
+yoff = 0.25*yspacing*ysign ; VERT. OFFSET FOR SYM/LINE.
+
+x0 = position[0] + (margin)*xspacing ; INITIAL X & Y POSITIONS
+y0 = position[1] - margin*yspacing + yalign*yspacing ; WELL, THIS WORKS!
+;
+; =====>> OUTPUT TEXT FOR LEGEND, ITEM BY ITEM.
+; =====>> FOR EACH ITEM, PLACE SYM/LINE, THEN DELIMITER,
+; =====>> THEN TEXT---UPDATING X & Y POSITIONS EACH TIME.
+; =====>> THERE ARE A NUMBER OF EXCEPTIONS DONE WITH IF STATEMENTS.
+;
+for iclr = 0,clear do begin
+ y = y0 ; STARTING X & Y POSITIONS
+ x = x0
+ if ltor then xend = 0 else xend = 1 ; SAVED WIDTH FOR DRAWING BOX
+
+ if ttob then ii = [0,n-1,1] else ii = [n-1,0,-1]
+ for i = ii[0],ii[1],ii[2] do begin
+ if vertical then x = x0 else y = y0 ; RESET EITHER X OR Y
+ x = x + xspacing ; UPDATE X & Y POSITIONS
+ y = y - yspacing
+ if nlpv eq 0 then goto,TEXT_ONLY ; FLAG FOR TEXT ONLY
+ if (psym[i] eq 0) and (vectorfont[i] eq '') then num = (number + 1) > 3 else num = number
+ if psym[i] lt 0 then num = number > 2 ; TO SHOW CONNECTING LINE
+ if psym[i] eq 0 then expand = 1 else expand = 2
+ xp = x + expand*pspacing*indgen(num)*xspacing
+ if (psym[i] gt 0) and (num eq 1) and vertical then xp = x + xt/2.
+ yp = y + intarr(num)
+ if vectorfont[i] eq '' then yp = yp + yoff
+ if psym[i] eq 0 then begin
+ xp = [min(xp),max(xp)] ; TO EXPOSE LINESTYLES
+ yp = [min(yp),max(yp)] ; DITTO
+ endif
+ if (psym[i] eq 8) and (N_elements(usersym) GT 1) then $
+ usersym,usersym,fill=fill,color=colors[i]
+;; extra by djseed .. psym=88 means use the already defined usersymbol
+ if psym[i] eq 88 then psym[i] =8
+ if vectorfont[i] ne '' then begin
+; if (num eq 1) and vertical then xp = x + xt/2 ; IF 1, CENTERED.
+ xyouts,xp,yp,vectorfont[i],width=width,color=colors[i] $
+ ,size=charsize,align=xalign,charthick = charthick,/norm
+ xt = xt > width
+ xp = xp + width/2.
+ endif else begin
+ if symline and (linestyle[i] ge 0) then plots,xp,yp,color=colors[i] $
+ ,/normal,linestyle=linestyle[i],psym=psym[i],symsize=symsize[i], $
+ thick=thick[i]
+ endelse
+
+ if vertical then x = x + xt else if ltor then x = max(xp) else x = min(xp)
+ if symline then x = x + xspacing
+ TEXT_ONLY:
+ if vertical and (vectorfont[i] eq '') and symline and (linestyle[i] eq -99) then x=x0 + xspacing
+ xyouts,x,y,delimiter,width=width,/norm,color=textcolors[i], $
+ size=charsize,align=xalign,charthick = charthick
+ x = x + width*xsign
+ if width ne 0 then x = x + 0.5*xspacing
+ xyouts,x,y,items[i],width=width,/norm,color=textcolors[i],size=charsize, $
+ align=xalign,charthick=charthick
+ x = x + width*xsign
+ if not vertical and (i lt (n-1)) then x = x+2*xspacing; ADD INTER-ITEM SPACE
+ xfinal = (x + xspacing*margin)
+ if ltor then xend = xfinal > xend else xend = xfinal < xend ; UPDATE END X
+ endfor
+
+ if (iclr lt clear ) then begin
+; =====>> CLEAR AREA
+ x = position[0]
+ y = position[1]
+ if vertical then bottom = n else bottom = 1
+ ywidth = - (2*margin+bottom-0.5)*yspacing
+ corners = [x,y+ywidth,xend,y]
+ polyfill,[x,xend,xend,x,x],y + [0,0,ywidth,ywidth,0],/norm,color=-1
+; plots,[x,xend,xend,x,x],y + [0,0,ywidth,ywidth,0],thick=2
+ endif else begin
+
+;
+; =====>> OUTPUT BORDER
+;
+ x = position[0]
+ y = position[1]
+ if vertical then bottom = n else bottom = 1
+ ywidth = - (2*margin+bottom-0.5)*yspacing
+ corners = [x,y+ywidth,xend,y]
+ if box then plots,[x,xend,xend,x,x],y + [0,0,ywidth,ywidth,0],/norm
+ return
+ endelse
+endfor
+
+end
+
diff --git a/examples/IsothermalPotential/GravityOnly/makeIC.py b/examples/IsothermalPotential/GravityOnly/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..88115668bd11e19f93765860540dcf33c6ae5c64
--- /dev/null
+++ b/examples/IsothermalPotential/GravityOnly/makeIC.py
@@ -0,0 +1,168 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ # Tom Theuns (tom.theuns@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+import h5py
+import sys
+import numpy
+import math
+import random
+
+# Generates N particles in a spherical distribution centred on [0,0,0], to be moved in an isothermal potential
+# usage: python makeIC.py 1000 0 : generate 1000 particles on circular orbits
+# python makeIC.py 1000 1 : generate 1000 particles with Lz/L uniform in [0,1]
+# all particles move in the xy plane, and start at y=0
+
+# physical constants in cgs
+NEWTON_GRAVITY_CGS = 6.672e-8
+SOLAR_MASS_IN_CGS = 1.9885e33
+PARSEC_IN_CGS = 3.0856776e18
+PROTON_MASS_IN_CGS = 1.6726231e24
+YEAR_IN_CGS = 3.154e+7
+
+# choice of units
+const_unit_length_in_cgs = (1000*PARSEC_IN_CGS)
+const_unit_mass_in_cgs = (SOLAR_MASS_IN_CGS)
+const_unit_velocity_in_cgs = (1e5)
+
+print "UnitMass_in_cgs: ", const_unit_mass_in_cgs
+print "UnitLength_in_cgs: ", const_unit_length_in_cgs
+print "UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs
+
+
+# rotation speed of isothermal potential [km/s]
+vrot_kms = 200.
+
+# derived units
+const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs)
+const_G = ((NEWTON_GRAVITY_CGS*const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/(const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs)))
+print 'G=', const_G
+vrot = vrot_kms * 1e5 / const_unit_velocity_in_cgs
+
+
+# Parameters
+periodic= 1 # 1 For periodic box
+boxSize = 400. # [kpc]
+Radius = 100. # maximum radius of particles [kpc]
+G = const_G
+
+N = int(sys.argv[1]) # Number of particles
+icirc = int(sys.argv[2]) # if = 0, all particles are on circular orbits, if = 1, Lz/Lcirc uniform in ]0,1[
+L = N**(1./3.)
+
+# these are not used but necessary for I/O
+rho = 2. # Density
+P = 1. # Pressure
+gamma = 5./3. # Gas adiabatic index
+fileName = "Isothermal.hdf5"
+
+
+#---------------------------------------------------
+numPart = N
+mass = 1
+internalEnergy = P / ((gamma - 1.)*rho)
+
+#--------------------------------------------------
+
+#File
+file = h5py.File(fileName, 'w')
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = const_unit_length_in_cgs
+grp.attrs["Unit mass in cgs (U_M)"] = const_unit_mass_in_cgs
+grp.attrs["Unit time in cgs (U_t)"] = const_unit_length_in_cgs / const_unit_velocity_in_cgs
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+grp.attrs["NumPart_Total"] = [0, numPart, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [0, numPart, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
+
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+# set seed for random number
+numpy.random.seed(1234)
+
+#Particle group
+#grp0 = file.create_group("/PartType0")
+grp1 = file.create_group("/PartType1")
+#generate particle positions
+radius = Radius * (numpy.random.rand(N))**(1./3.)
+ctheta = -1. + 2 * numpy.random.rand(N)
+stheta = numpy.sqrt(1.-ctheta**2)
+phi = 2 * math.pi * numpy.random.rand(N)
+r = numpy.zeros((numPart, 3))
+#r[:,0] = radius * stheta * numpy.cos(phi)
+#r[:,1] = radius * stheta * numpy.sin(phi)
+#r[:,2] = radius * ctheta
+r[:,0] = radius
+#
+speed = vrot
+v = numpy.zeros((numPart, 3))
+omega = speed / radius
+period = 2.*math.pi/omega
+print 'period = minimum = ',min(period), ' maximum = ',max(period)
+
+omegav = omega
+if (icirc != 0):
+ omegav = omega * numpy.random.rand(N)
+
+v[:,0] = -omegav * r[:,1]
+v[:,1] = omegav * r[:,0]
+
+ds = grp1.create_dataset('Velocities', (numPart, 3), 'f')
+ds[()] = v
+v = numpy.zeros(1)
+
+m = numpy.full((numPart, ), mass)
+ds = grp1.create_dataset('Masses', (numPart,), 'f')
+ds[()] = m
+m = numpy.zeros(1)
+
+h = numpy.full((numPart, ), 1.1255 * boxSize / L)
+ds = grp1.create_dataset('SmoothingLength', (numPart,), 'f')
+ds[()] = h
+h = numpy.zeros(1)
+
+u = numpy.full((numPart, ), internalEnergy)
+ds = grp1.create_dataset('InternalEnergy', (numPart,), 'f')
+ds[()] = u
+u = numpy.zeros(1)
+
+
+ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False, dtype='L')
+ds = grp1.create_dataset('ParticleIDs', (numPart, ), 'L')
+ds[()] = ids
+
+ds = grp1.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = r
+
+
+file.close()
diff --git a/examples/IsothermalPotential/GravityOnly/readme.txt b/examples/IsothermalPotential/GravityOnly/readme.txt
new file mode 100644
index 0000000000000000000000000000000000000000..90fb1872aa2301cab133b8a20a8cd8de724d4553
--- /dev/null
+++ b/examples/IsothermalPotential/GravityOnly/readme.txt
@@ -0,0 +1,6 @@
+;
+; this probelm generates a set of gravity particles in an isothermal
+; potential and follows their orbits. Tests verify consdevation of
+; energy and angular momentum
+;
+;
diff --git a/examples/IsothermalPotential/GravityOnly/run.sh b/examples/IsothermalPotential/GravityOnly/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..89e7c10ce19ce0576127551f63eea5c9e4e955ff
--- /dev/null
+++ b/examples/IsothermalPotential/GravityOnly/run.sh
@@ -0,0 +1,10 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e Isothermal.hdf5 ]
+then
+ echo "Generating initial conditions for the point mass potential box example..."
+ python makeIC.py 1000
+fi
+
+../swift -g -t 2 externalGravity.yml
diff --git a/examples/IsothermalPotential/GravityOnly/test.pro b/examples/IsothermalPotential/GravityOnly/test.pro
new file mode 100644
index 0000000000000000000000000000000000000000..ac1f5e915c38c107408ffed9579c52944295f079
--- /dev/null
+++ b/examples/IsothermalPotential/GravityOnly/test.pro
@@ -0,0 +1,168 @@
+;
+; test energy / angular momentum conservation of test problem
+;
+
+iplot = 1 ; if iplot = 1, make plot of E/Lz conservation, else, simply compare final and initial energy
+
+; set physical constants
+@physunits
+
+indir = './'
+basefile = 'Isothermal_'
+
+; set properties of potential
+uL = 1e3 * phys.pc ; unit of length
+uM = phys.msun ; unit of mass
+uV = 1d5 ; unit of velocity
+vrot = 200. ; km/s
+r200 = 100. ; virial radius
+
+; derived units
+constG = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ;
+pcentre = [100.,100.,100.] * 1d3 * pc / uL
+
+;
+infile = indir + basefile + '*'
+spawn,'ls -1 '+infile,res
+nfiles = n_elements(res)
+
+
+
+; choose: calculate change of energy and Lz, comparing first and last
+; snapshots for all particles, or do so for a subset
+
+; compare all
+ifile = 0
+inf = indir + basefile + strtrim(string(ifile,'(i3.3)'),1) + '.hdf5'
+id = h5rd(inf,'PartType1/ParticleIDs')
+nfollow = n_elements(id)
+
+; follow a subset
+nfollow = 500 ; number of particles to follow
+
+;
+if (iplot eq 1) then begin
+ nskip = 1
+ nsave = nfiles
+endif else begin
+ nskip = nfiles - 2
+ nsave = 2
+endelse
+
+;
+lout = fltarr(nfollow, nsave) ; Lz
+xout = fltarr(nfollow, nsave) ; x
+yout = fltarr(nfollow, nsave) ; y
+zout = fltarr(nfollow, nsave) ; z
+eout = fltarr(nfollow, nsave) ; energies
+ekin = fltarr(nfollow, nsave)
+epot = fltarr(nfollow, nsave)
+tout = fltarr(nsave)
+
+
+
+ifile = 0
+isave = 0
+for ifile=0,nfiles-1,nskip do begin
+ inf = indir + basefile + strtrim(string(ifile,'(i3.3)'),1) + '.hdf5'
+ time = h5ra(inf, 'Header','Time')
+ p = h5rd(inf,'PartType1/Coordinates')
+ v = h5rd(inf,'PartType1/Velocities')
+ id = h5rd(inf,'PartType1/ParticleIDs')
+ indx = sort(id)
+;
+ id = id[indx]
+ for ic=0,2 do begin
+ tmp = reform(p[ic,*]) & p[ic,*] = tmp[indx]
+ tmp = reform(v[ic,*]) & v[ic,*] = tmp[indx]
+ endfor
+
+
+; calculate energy
+ dd = size(p,/dimen) & npart = dd[1]
+ ener = fltarr(npart)
+ dr = fltarr(npart) & dv = dr
+ for ic=0,2 do dr[*] = dr[*] + (p[ic,*]-pcentre[ic])^2
+ for ic=0,2 do dv[*] = dv[*] + v[ic,*]^2
+ xout[*,isave] = p[0,0:nfollow-1]-pcentre[0]
+ yout[*,isave] = p[1,0:nfollow-1]-pcentre[1]
+ zout[*,isave] = p[2,0:nfollow-1]-pcentre[2]
+ Lz = (p[0,*]-pcentre[0]) * v[1,*] - (p[1,*]-pcentre[1]) * v[0,*]
+ dr = sqrt(dr)
+; print,'time = ',time,p[0,0],v[0,0],id[0]
+ ek = 0.5 * dv
+; ep = - constG * mextern / dr
+ ep = -vrot*vrot * (1 + alog(r200/dr))
+ ener = ek + ep
+ tout(isave) = time
+ lout[*,isave] = lz[0:nfollow-1]
+ eout(*,isave) = ener[0:nfollow-1]
+ ekin(*,isave) = ek[0:nfollow-1]
+ epot(*,isave) = ep[0:nfollow-1]
+
+; write some output
+; print,' time= ',time,' e= ',eout[0],' Lz= ',lz[0],format='(%a %f %a
+; %f)'
+ print,format='('' time= '',f7.1,'' E= '',f9.2,'' Lz= '',e9.2)', time,eout[0],lz[0]
+ isave = isave + 1
+
+endfor
+x0 = reform(xout[0,*])
+y0 = reform(xout[1,*])
+z0 = reform(xout[2,*])
+
+; calculate relative energy change
+de = 0.0 * eout
+dl = 0.0 * lout
+nsave = isave
+for ifile=1, nsave-1 do de[*,ifile] = (eout[*,ifile]-eout[*,0])/eout[*,0]
+for ifile=1, nsave-1 do dl[*,ifile] = (lout[*,ifile] - lout[*,0])/lout[*,0]
+
+
+; calculate statistics of energy changes
+print,' relatve energy change: (per cent) ',minmax(de) * 100.
+print,' relative Lz change: (per cent) ',minmax(dl) * 100.
+
+; plot enery and Lz conservation for some particles
+if(iplot eq 1) then begin
+; plot results on energy conservation for some particles
+ nplot = min(10, nfollow)
+ wset,0
+ xr = [min(tout), max(tout)]
+ yr = [-2,2]*1d-2 ; in percent
+ plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/nodata,xtitle='time',ytitle='dE/E, dL/L (%)'
+ for i=0,nplot-1 do oplot,tout,de[i,*]
+ for i=0,nplot-1 do oplot,tout,dl[i,*],color=red
+ legend,['dE/E','dL/L'],linestyle=[0,0],color=[black,red],box=0,/bottom,/left
+ screen_to_png,'e-time.png'
+
+; plot orbits of those particles
+ wset,2
+ xr = [-100,100]
+ yr = xr
+ plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/iso,/nodata,xtitle='x',ytitle='y'
+ color = floor(findgen(nplot)*255/float(nplot))
+ for i=0,nplot-1 do oplot,xout[i,*],yout[i,*],color=color(i)
+ screen_to_png,'orbit.png'
+
+; plot radial position of these particles
+ wset,4
+ xr = [min(tout), max(tout)]
+ yr = [0,80]
+ plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/nodata,xtitle='t',ytitle='r'
+ color = floor(findgen(nplot)*255/float(nplot))
+for i=0,nplot-1 do begin dr = sqrt(reform(xout[i,*])^2 + reform(yout[i,*])^2) & oplot,tout,dr,color=color[i] & endfor
+ screen_to_png,'r-time.png'
+
+; make histogram of energy changes at end
+ wset,6
+ ohist,de,x,y,-0.05,0.05,0.001
+ plot,x,y,psym=10,xtitle='de (%)'
+ screen_to_png,'de-hist.png'
+
+
+endif
+
+end
+
+
diff --git a/src/const.h b/src/const.h
index 65404d6a3e89d85c64c81e28c6eaf9e1884c82e2..0487da30ccbd4fb7951de62a388c917b1d5510f3 100644
--- a/src/const.h
+++ b/src/const.h
@@ -60,7 +60,8 @@
//#define DEFAULT_SPH
/* Gravity properties */
-#define EXTERNAL_POTENTIAL_POINTMASS
+/* #define EXTERNAL_POTENTIAL_POINTMASS */
+#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
/* Are we debugging ? */
//#define SWIFT_DEBUG_CHECKS
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 7ef6db76e47f5a1cb2efee9d90fabb9f2ca7dd7f..0e4784693ad05b64e876b0b3f0f87b842a6efe9f 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -42,6 +42,9 @@ gravity_compute_timestep_external(const struct external_potential* potential,
dt =
fminf(dt, external_gravity_pointmass_timestep(potential, phys_const, gp));
#endif
+#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+ dt = fminf(dt, external_gravity_isothermalpotential_timestep(potential, phys_const, gp));
+#endif
return dt;
}
@@ -117,6 +120,9 @@ __attribute__((always_inline)) INLINE static void external_gravity(
#ifdef EXTERNAL_POTENTIAL_POINTMASS
external_gravity_pointmass(potential, phys_const, gp);
#endif
+#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+ external_gravity_isothermalpotential(potential, phys_const, gp);
+#endif
}
/**
diff --git a/src/potentials.c b/src/potentials.c
index 98e57d7959963e07cfbbd2fdb60a0df504e6a487..d55d3d669f1a146a0932c46b551ee06f37ed3921 100644
--- a/src/potentials.c
+++ b/src/potentials.c
@@ -48,6 +48,16 @@ void potential_init(const struct swift_params* parameter_file,
parser_get_param_double(parameter_file, "PointMass:mass");
#endif /* EXTERNAL_POTENTIAL_POINTMASS */
+#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+ potential -> isothermal_potential.x =
+ parser_get_param_double(parameter_file, "IsothermalPotential:position_x");
+ potential -> isothermal_potential.y =
+ parser_get_param_double(parameter_file, "IsothermalPotential:position_y");
+ potential -> isothermal_potential.z =
+ parser_get_param_double(parameter_file, "IsothermalPotential:position_z");
+ potential -> isothermal_potential.vrot =
+ parser_get_param_double(parameter_file, "IsothermalPotential:vrot");
+#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
}
/**
@@ -62,4 +72,9 @@ void potential_print(const struct external_potential* potential) {
potential->point_mass.x, potential->point_mass.y,
potential->point_mass.z, potential->point_mass.mass);
#endif /* EXTERNAL_POTENTIAL_POINTMASS */
+#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+ message("Isothermal potential properties are (x,y,z) = (%e, %e, %e), vrot = %e",
+ potential->isothermal_potential.x, potential->isothermal_potential.y,
+ potential->isothermal_potential.z, potential->isothermal_potential.vrot);
+#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
}
diff --git a/src/potentials.h b/src/potentials.h
index c3db5a3d2231c93b1d83aa94bca83f7f7d7106bc..c2ee511af3e390a4942c176867f64b8c7bd43fa5 100644
--- a/src/potentials.h
+++ b/src/potentials.h
@@ -44,7 +44,68 @@ struct external_potential {
double mass;
} point_mass;
#endif
+#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+ struct {
+ double x, y, z;
+ double vrot;
+ } isothermal_potential;
+#endif
};
+/* Include external isothermal potential */
+#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+
+#define EXTERNAL_GRAVITY_TIMESTEP_PREFACTOR 0.03f
+/**
+ * @brief Computes the time-step due to the acceleration from a point mass
+ *
+ * @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(
+ const struct external_potential* potential,
+ const struct phys_const* const phys_const,
+ const struct gpart* const g) {
+
+ const float dx = g->x[0] - potential->isothermal_potential.x;
+ const float dy = g->x[1] - potential->isothermal_potential.y;
+ const float dz = g->x[2] - potential->isothermal_potential.z;
+ const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
+ const float drdv = dx * (g->v_full[0]) + dy * (g->v_full[1]) + dz * (g->v_full[2]);
+ const double vrot = potential->isothermal_potential.vrot;
+
+ const float dota_x = vrot * vrot * rinv2 * (g->v_full[0] - 2 * drdv * dx * rinv2);
+ const float dota_y = vrot * vrot * rinv2 * (g->v_full[1] - 2 * drdv * dy * rinv2);
+ const float dota_z = vrot * vrot * rinv2 * (g->v_full[2] - 2 * drdv * dz * rinv2);
+ const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
+ const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] +
+ g->a_grav[2] * g->a_grav[2];
+
+ // error(" dx= %e dvx= %e rinv2= %e a2= %e dota_2= %e dt= %e", dx, g->v_full[1], rinv2, a_2, dota_2, 0.03f * sqrtf(a_2 / dota_2));
+
+ return EXTERNAL_GRAVITY_TIMESTEP_PREFACTOR * sqrtf(a_2 / dota_2);
+}
+/**
+ * @brief Computes the gravitational acceleration of a particle due to a point
+ * mass
+ *
+ * @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) {
+
+ const float dx = g->x[0] - potential->isothermal_potential.x;
+ const float dy = g->x[1] - potential->isothermal_potential.y;
+ const float dz = g->x[2] - potential->isothermal_potential.z;
+ const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
+
+ const double vrot = potential->isothermal_potential.vrot;
+ g->a_grav[0] += -vrot * vrot * rinv2 * dx;
+ g->a_grav[1] += -vrot * vrot * rinv2 * dy;
+ g->a_grav[2] += -vrot * vrot * rinv2 * dz;
+ // error(" %f %f %f %f", vrot, rinv2, dx, g->a_grav[0]);
+}
+#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
/* Include exteral pointmass potential */
#ifdef EXTERNAL_POTENTIAL_POINTMASS