diff --git a/examples/Disk-Patch/GravityOnly/readme.txt b/examples/Disk-Patch/GravityOnly/README
similarity index 63%
rename from examples/Disk-Patch/GravityOnly/readme.txt
rename to examples/Disk-Patch/GravityOnly/README
index a18365146aa13abb81e7d3b8fcac31b92959ae3f..00673d274a50f9ab3043adb35d54629e1a993f0f 100644
--- a/examples/Disk-Patch/GravityOnly/readme.txt
+++ b/examples/Disk-Patch/GravityOnly/README
@@ -1,5 +1,5 @@
-sets-up for a potential of a patch f disk, see Creasey, Theuns &
-Bower,2013.
+Setup for a potential of a patch disk, see Creasey, Theuns &
+Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
rho(z) = (Sigma/2b) / cosh^2(z/b)
where Sigma is the surface density, and b the scale height
diff --git a/examples/Disk-Patch/GravityOnly/legend.pro b/examples/Disk-Patch/GravityOnly/legend.pro
deleted file mode 100755
index e510bdfcd7dcabed796cb090c6fa7d54536608b6..0000000000000000000000000000000000000000
--- a/examples/Disk-Patch/GravityOnly/legend.pro
+++ /dev/null
@@ -1,459 +0,0 @@
-;+
-; NAME:
-; LEGEND
-; PURPOSE:
-; Create an annotation legend for a plot.
-; EXPLANATION:
-; This procedure makes a legend for a plot. The legend can contain
-; a mixture of symbols, linestyles, Hershey characters (vectorfont),
-; and filled polygons (usersym). A test procedure, legendtest.pro,
-; shows legend's capabilities. Placement of the legend is controlled
-; with keywords like /right, /top, and /center or by using a position
-; keyword for exact placement (position=[x,y]) or via mouse (/position).
-; CALLING SEQUENCE:
-; LEGEND [,items][,keyword options]
-; EXAMPLES:
-; The call:
-; legend,['Plus sign','Asterisk','Period'],psym=[1,2,3]
-; produces:
-; -----------------
-; | |
-; | + Plus sign |
-; | * Asterisk |
-; | . Period |
-; | |
-; -----------------
-; Each symbol is drawn with a plots command, so they look OK.
-; Other examples are given in optional output keywords.
-;
-; lines = indgen(6) ; for line styles
-; items = 'linestyle '+strtrim(lines,2) ; annotations
-; legend,items,linestyle=lines ; vertical legend---upper left
-; items = ['Plus sign','Asterisk','Period']
-; sym = [1,2,3]
-; legend,items,psym=sym ; ditto except using symbols
-; legend,items,psym=sym,/horizontal ; horizontal format
-; legend,items,psym=sym,box=0 ; sans border
-; legend,items,psym=sym,delimiter='=' ; embed '=' betw psym & text
-; legend,items,psym=sym,margin=2 ; 2-character margin
-; legend,items,psym=sym,position=[x,y] ; upper left in data coords
-; legend,items,psym=sym,pos=[x,y],/norm ; upper left in normal coords
-; legend,items,psym=sym,pos=[x,y],/device ; upper left in device coords
-; legend,items,psym=sym,/position ; interactive position
-; legend,items,psym=sym,/right ; at upper right
-; legend,items,psym=sym,/bottom ; at lower left
-; legend,items,psym=sym,/center ; approximately near center
-; legend,items,psym=sym,number=2 ; plot two symbols, not one
-; legend,items,/fill,psym=[8,8,8],colors=[10,20,30]; 3 filled squares
-; INPUTS:
-; items = text for the items in the legend, a string array.
-; For example, items = ['diamond','asterisk','square'].
-; You can omit items if you don't want any text labels.
-; OPTIONAL INPUT KEYWORDS:
-;
-; linestyle = array of linestyle numbers If linestyle[i] < 0, then omit
-; ith symbol or line to allow a multi-line entry. If
-; linestyle = -99 then text will be left-justified.
-; psym = array of plot symbol numbers. If psym[i] is negative, then a
-; line connects pts for ith item. If psym[i] = 8, then the
-; procedure usersym is called with vertices define in the
-; keyword usersym. If psym[i] = 88, then use the previously
-; defined user symbol
-; vectorfont = vector-drawn characters for the sym/line column, e.g.,
-; ['!9B!3','!9C!3','!9D!3'] produces an open square, a checkmark,
-; and a partial derivative, which might have accompanying items
-; ['BOX','CHECK','PARTIAL DERIVATIVE'].
-; There is no check that !p.font is set properly, e.g., -1 for
-; X and 0 for PostScript. This can produce an error, e.g., use
-; !20 with PostScript and !p.font=0, but allows use of Hershey
-; *AND* PostScript fonts together.
-; N. B.: Choose any of linestyle, psym, and/or vectorfont. If none is
-; present, only the text is output. If more than one
-; is present, all need the same number of elements, and normal
-; plot behaviour occurs.
-; By default, if psym is positive, you get one point so there is
-; no connecting line. If vectorfont[i] = '',
-; then plots is called to make a symbol or a line, but if
-; vectorfont[i] is a non-null string, then xyouts is called.
-; /help = flag to print header
-; /horizontal = flag to make the legend horizontal
-; /vertical = flag to make the legend vertical (D=vertical)
-; box = flag to include/omit box around the legend (D=include)
-; clear = flag to clear the box area before drawing the legend
-; delimiter = embedded character(s) between symbol and text (D=none)
-; colors = array of colors for plot symbols/lines (D=!P.color)
-; textcolors = array of colors for text (D=!P.color)
-; margin = margin around text measured in characters and lines
-; spacing = line spacing (D=bit more than character height)
-; pspacing = psym spacing (D=3 characters) (when number of symbols is
-; greater than 1)
-; charsize = just like !p.charsize for plot labels
-; charthick = just like !p.charthick for plot labels
-; thick = array of line thickness numbers (D = !P.thick), if used, then
-; linestyle must also be specified
-; position = data coordinates of the /top (D) /left (D) of the legend
-; normal = use normal coordinates for position, not data
-; device = use device coordinates for position, not data
-; number = number of plot symbols to plot or length of line (D=1)
-; usersym = 2-D array of vertices, cf. usersym in IDL manual.
-; (/USERSYM =square, default is to use existing USERSYM definition)
-; /fill = flag to fill the usersym
-; /left_legend = flag to place legend snug against left side of plot
-; window (D)
-; /right_legend = flag to place legend snug against right side of plot
-; window. If /right,pos=[x,y], then x is position of RHS and
-; text runs right-to-left.
-; /top_legend = flag to place legend snug against top of plot window (D)
-; /bottom = flag to place legend snug against bottom of plot window
-; /top,pos=[x,y] and /bottom,pos=[x,y] produce same positions.
-;
-; If LINESTYLE, PSYM, VECTORFONT, THICK, COLORS, or TEXTCOLORS are
-; supplied as scalars, then the scalar value is set for every line or
-; symbol in the legend.
-; Outputs:
-; legend to current plot device
-; OPTIONAL OUTPUT KEYWORDS:
-; corners = 4-element array, like !p.position, of the normalized
-; coords for the box (even if box=0): [llx,lly,urx,ury].
-; Useful for multi-column or multi-line legends, for example,
-; to make a 2-column legend, you might do the following:
-; c1_items = ['diamond','asterisk','square']
-; c1_psym = [4,2,6]
-; c2_items = ['solid','dashed','dotted']
-; c2_line = [0,2,1]
-; legend,c1_items,psym=c1_psym,corners=c1,box=0
-; legend,c2_items,line=c2_line,corners=c2,box=0,pos=[c1[2],c1[3]]
-; c = [c1[0]<c2[0],c1[1]<c2[1],c1[2]>c2[2],c1[3]>c2[3]]
-; plots,[c[0],c[0],c[2],c[2],c[0]],[c[1],c[3],c[3],c[1],c[1]],/norm
-; Useful also to place the legend. Here's an automatic way to place
-; the legend in the lower right corner. The difficulty is that the
-; legend's width is unknown until it is plotted. In this example,
-; the legend is plotted twice: the first time in the upper left, the
-; second time in the lower right.
-; legend,['1','22','333','4444'],linestyle=indgen(4),corners=corners
-; ; BOGUS LEGEND---FIRST TIME TO REPORT CORNERS
-; xydims = [corners[2]-corners[0],corners[3]-corners[1]]
-; ; SAVE WIDTH AND HEIGHT
-; chdim=[!d.x_ch_size/float(!d.x_size),!d.y_ch_size/float(!d.y_size)]
-; ; DIMENSIONS OF ONE CHARACTER IN NORMALIZED COORDS
-; pos = [!x.window[1]-chdim[0]-xydims[0] $
-; ,!y.window[0]+chdim[1]+xydims[1]]
-; ; CALCULATE POSITION FOR LOWER RIGHT
-; plot,findgen(10) ; SIMPLE PLOT; YOU DO WHATEVER YOU WANT HERE.
-; legend,['1','22','333','4444'],linestyle=indgen(4),pos=pos
-; ; REDO THE LEGEND IN LOWER RIGHT CORNER
-; You can modify the pos calculation to place the legend where you
-; want. For example to place it in the upper right:
-; pos = [!x.window[1]-chdim[0]-xydims[0],!y.window[1]-xydims[1]]
-; Common blocks:
-; none
-; Procedure:
-; If keyword help is set, call doc_library to print header.
-; See notes in the code. Much of the code deals with placement of the
-; legend. The main problem with placement is not being
-; able to sense the length of a string before it is output. Some crude
-; approximations are used for centering.
-; Restrictions:
-; Here are some things that aren't implemented.
-; - An orientation keyword would allow lines at angles in the legend.
-; - An array of usersyms would be nice---simple change.
-; - An order option to interchange symbols and text might be nice.
-; - Somebody might like double boxes, e.g., with box = 2.
-; - Another feature might be a continuous bar with ticks and text.
-; - There are no guards to avoid writing outside the plot area.
-; - There is no provision for multi-line text, e.g., '1st line!c2nd line'
-; Sensing !c would be easy, but !c isn't implemented for PostScript.
-; A better way might be to simply output the 2nd line as another item
-; but without any accompanying symbol or linestyle. A flag to omit
-; the symbol and linestyle is linestyle[i] = -1.
-; - There is no ability to make a title line containing any of titles
-; for the legend, for the symbols, or for the text.
-; Side Effects:
-; Modification history:
-; write, 24-25 Aug 92, F K Knight (knight@ll.mit.edu)
-; allow omission of items or omission of both psym and linestyle, add
-; corners keyword to facilitate multi-column legends, improve place-
-; ment of symbols and text, add guards for unequal size, 26 Aug 92, FKK
-; add linestyle(i)=-1 to suppress a single symbol/line, 27 Aug 92, FKK
-; add keyword vectorfont to allow characters in the sym/line column,
-; 28 Aug 92, FKK
-; add /top, /bottom, /left, /right keywords for automatic placement at
-; the four corners of the plot window. The /right keyword forces
-; right-to-left printing of menu. 18 Jun 93, FKK
-; change default position to data coords and add normal, data, and
-; device keywords, 17 Jan 94, FKK
-; add /center keyword for positioning, but it is not precise because
-; text string lengths cannot be known in advance, 17 Jan 94, FKK
-; add interactive positioning with /position keyword, 17 Jan 94, FKK
-; allow a legend with just text, no plotting symbols. This helps in
-; simply describing a plot or writing assumptions done, 4 Feb 94, FKK
-; added thick, symsize, and clear keyword Feb 96, W. Landsman HSTX
-; David Seed, HR Wallingford, d.seed@hrwallingford.co.uk
-; allow scalar specification of keywords, Mar 96, W. Landsman HSTX
-; added charthick keyword, June 96, W. Landsman HSTX
-; Made keyword names left,right,top,bottom,center longer,
-; Aug 16, 2000, Kim Tolbert
-; Added ability to have regular text lines in addition to plot legend
-; lines in legend. If linestyle is -99 that item is left-justified.
-; Previously, only option for no sym/line was linestyle=-1, but then text
-; was lined up after sym/line column. 10 Oct 2000, Kim Tolbert
-; Make default value of thick = !P.thick W. Landsman Jan. 2001
-; Don't overwrite existing USERSYM definition W. Landsman Mar. 2002
-;-
-pro legend, items, BOTTOM_LEGEND=bottom, BOX = box, CENTER_LEGEND=center, $
- CHARTHICK=charthick, CHARSIZE = charsize, CLEAR = clear, COLORS = colorsi, $
- CORNERS = corners, DATA=data, DELIMITER=delimiter, DEVICE=device, $
- FILL=fill, HELP = help, HORIZONTAL=horizontal,LEFT_LEGEND=left, $
- LINESTYLE=linestylei, MARGIN=margin, NORMAL=normal, NUMBER=number, $
- POSITION=position,PSPACING=pspacing, PSYM=psymi, RIGHT_LEGEND=right, $
- SPACING=spacing, SYMSIZE=symsize, TEXTCOLORS=textcolorsi, THICK=thicki, $
- TOP_LEGEND=top, USERSYM=usersym, VECTORFONT=vectorfonti, VERTICAL=vertical
-;
-; =====>> HELP
-;
-on_error,2
-if keyword_set(help) then begin & doc_library,'legend' & return & endif
-;
-; =====>> SET DEFAULTS FOR SYMBOLS, LINESTYLES, AND ITEMS.
-;
- ni = n_elements(items)
- np = n_elements(psymi)
- nl = n_elements(linestylei)
- nth = n_elements(thicki)
- nv = n_elements(vectorfonti)
- nlpv = max([np,nl,nv])
- n = max([ni,np,nl,nv]) ; NUMBER OF ENTRIES
-strn = strtrim(n,2) ; FOR ERROR MESSAGES
-if n eq 0 then message,'No inputs! For help, type legend,/help.'
-if ni eq 0 then begin
- items = replicate('',n) ; DEFAULT BLANK ARRAY
-endif else begin
- if size(items,/TNAME) NE 'STRING' then message, $
- 'First parameter must be a string array. For help, type legend,/help.'
- if ni ne n then message,'Must have number of items equal to '+strn
-endelse
-symline = (np ne 0) or (nl ne 0) ; FLAG TO PLOT SYM/LINE
- if (np ne 0) and (np ne n) and (np NE 1) then message, $
- 'Must have 0, 1 or '+strn+' elements in PSYM array.'
- if (nl ne 0) and (nl ne n) and (nl NE 1) then message, $
- 'Must have 0, 1 or '+strn+' elements in LINESTYLE array.'
- if (nth ne 0) and (nth ne n) and (nth NE 1) then message, $
- 'Must have 0, 1 or '+strn+' elements in THICK array.'
-
- case nl of
- 0: linestyle = intarr(n) ;Default = solid
- 1: linestyle = intarr(n) + linestylei
- else: linestyle = linestylei
- endcase
-
- case nth of
- 0: thick = replicate(!p.thick,n) ;Default = !P.THICK
- 1: thick = intarr(n) + thicki
- else: thick = thicki
- endcase
-
- case np of ;Get symbols
- 0: psym = intarr(n) ;Default = solid
- 1: psym = intarr(n) + psymi
- else: psym = psymi
- endcase
-
- case nv of
- 0: vectorfont = replicate('',n)
- 1: vectorfont = replicate(vectorfonti,n)
- else: vectorfont = vectorfonti
- endcase
-;
-; =====>> CHOOSE VERTICAL OR HORIZONTAL ORIENTATION.
-;
-if n_elements(horizontal) eq 0 then begin ; D=VERTICAL
- if n_elements(vertical) eq 0 then vertical = 1
-endif else begin
- if n_elements(vertical) eq 0 then vertical = not horizontal
-endelse
-;
-; =====>> SET DEFAULTS FOR OTHER OPTIONS.
-;
-if n_elements(box) eq 0 then box = 1
-if n_elements(clear) eq 0 then clear = 0
-
-if n_elements(margin) eq 0 then margin = 0.5
-if n_elements(delimiter) eq 0 then delimiter = ''
-if n_elements(charsize) eq 0 then charsize = !p.charsize
-if n_elements(charthick) eq 0 then charthick = !p.charthick
-if charsize eq 0 then charsize = 1
-if (n_elements (symsize) eq 0) then symsize= charsize + intarr(n)
-if n_elements(number) eq 0 then number = 1
- case N_elements(colorsi) of
- 0: colors = replicate(!P.color,n) ;Default is !P.COLOR
- 1: colors = replicate(colorsi,n)
- else: colors = colorsi
- endcase
-
- case N_elements(textcolorsi) of
- 0: textcolors = replicate(!P.color,n) ;Default is !P.COLOR
- 1: textcolors = replicate(textcolorsi,n)
- else: textcolors = textcolorsi
- endcase
- fill = keyword_set(fill)
-if n_elements(usersym) eq 1 then usersym = 2*[[0,0],[0,1],[1,1],[1,0],[0,0]]-1
-;
-; =====>> INITIALIZE SPACING
-;
-if n_elements(spacing) eq 0 then spacing = 1.2
-if n_elements(pspacing) eq 0 then pspacing = 3
-xspacing = !d.x_ch_size/float(!d.x_size) * (spacing > charsize)
-yspacing = !d.y_ch_size/float(!d.y_size) * (spacing > charsize)
-ltor = 1 ; flag for left-to-right
-if n_elements(left) eq 1 then ltor = left eq 1
-if n_elements(right) eq 1 then ltor = right ne 1
-ttob = 1 ; flag for top-to-bottom
-if n_elements(top) eq 1 then ttob = top eq 1
-if n_elements(bottom) eq 1 then ttob = bottom ne 1
-xalign = ltor ne 1 ; x alignment: 1 or 0
-yalign = -0.5*ttob + 1 ; y alignment: 0.5 or 1
-xsign = 2*ltor - 1 ; xspacing direction: 1 or -1
-ysign = 2*ttob - 1 ; yspacing direction: 1 or -1
-if not ttob then yspacing = -yspacing
-if not ltor then xspacing = -xspacing
-;
-; =====>> INITIALIZE POSITIONS: FIRST CALCULATE X OFFSET FOR TEXT
-;
-xt = 0
-if nlpv gt 0 then begin ; SKIP IF TEXT ITEMS ONLY.
-if vertical then begin ; CALC OFFSET FOR TEXT START
- for i = 0,n-1 do begin
- if (psym[i] eq 0) and (vectorfont[i] eq '') then num = (number + 1) > 3 else num = number
- if psym[i] lt 0 then num = number > 2 ; TO SHOW CONNECTING LINE
- if psym[i] eq 0 then expand = 1 else expand = 2
- thisxt = (expand*pspacing*(num-1)*xspacing)
- if ltor then xt = thisxt > xt else xt = thisxt < xt
- endfor
-endif ; NOW xt IS AN X OFFSET TO ALIGN ALL TEXT ENTRIES.
-endif
-;
-; =====>> INITIALIZE POSITIONS: SECOND LOCATE BORDER
-;
-if !x.window[0] eq !x.window[1] then begin
- plot,/nodata,xstyle=4,ystyle=4,[0],/noerase
-endif
-; next line takes care of weirdness with small windows
-pos = [min(!x.window),min(!y.window),max(!x.window),max(!y.window)]
-case n_elements(position) of
- 0: begin
- if ltor then px = pos[0] else px = pos[2]
- if ttob then py = pos[3] else py = pos[1]
- if keyword_set(center) then begin
- if not keyword_set(right) and not keyword_set(left) then $
- px = (pos[0] + pos[2])/2. - xt
- if not keyword_set(top) and not keyword_set(bottom) then $
- py = (pos[1] + pos[3])/2. + n*yspacing
- endif
- position = [px,py] + [xspacing,-yspacing]
- end
- 1: begin ; interactive
- message,/inform,'Place mouse at upper left corner and click any mouse button.'
- cursor,x,y,/normal
- position = [x,y]
- end
- 2: begin ; convert upper left corner to normal coordinates
- if keyword_set(data) then $
- position = convert_coord(position,/to_norm) $
- else if keyword_set(device) then $
- position = convert_coord(position,/to_norm,/device) $
- else if not keyword_set(normal) then $
- position = convert_coord(position,/to_norm)
- end
- else: message,'Position keyword can have 0, 1, or 2 elements only. Try legend,/help.'
-endcase
-
-yoff = 0.25*yspacing*ysign ; VERT. OFFSET FOR SYM/LINE.
-
-x0 = position[0] + (margin)*xspacing ; INITIAL X & Y POSITIONS
-y0 = position[1] - margin*yspacing + yalign*yspacing ; WELL, THIS WORKS!
-;
-; =====>> OUTPUT TEXT FOR LEGEND, ITEM BY ITEM.
-; =====>> FOR EACH ITEM, PLACE SYM/LINE, THEN DELIMITER,
-; =====>> THEN TEXT---UPDATING X & Y POSITIONS EACH TIME.
-; =====>> THERE ARE A NUMBER OF EXCEPTIONS DONE WITH IF STATEMENTS.
-;
-for iclr = 0,clear do begin
- y = y0 ; STARTING X & Y POSITIONS
- x = x0
- if ltor then xend = 0 else xend = 1 ; SAVED WIDTH FOR DRAWING BOX
-
- if ttob then ii = [0,n-1,1] else ii = [n-1,0,-1]
- for i = ii[0],ii[1],ii[2] do begin
- if vertical then x = x0 else y = y0 ; RESET EITHER X OR Y
- x = x + xspacing ; UPDATE X & Y POSITIONS
- y = y - yspacing
- if nlpv eq 0 then goto,TEXT_ONLY ; FLAG FOR TEXT ONLY
- if (psym[i] eq 0) and (vectorfont[i] eq '') then num = (number + 1) > 3 else num = number
- if psym[i] lt 0 then num = number > 2 ; TO SHOW CONNECTING LINE
- if psym[i] eq 0 then expand = 1 else expand = 2
- xp = x + expand*pspacing*indgen(num)*xspacing
- if (psym[i] gt 0) and (num eq 1) and vertical then xp = x + xt/2.
- yp = y + intarr(num)
- if vectorfont[i] eq '' then yp = yp + yoff
- if psym[i] eq 0 then begin
- xp = [min(xp),max(xp)] ; TO EXPOSE LINESTYLES
- yp = [min(yp),max(yp)] ; DITTO
- endif
- if (psym[i] eq 8) and (N_elements(usersym) GT 1) then $
- usersym,usersym,fill=fill,color=colors[i]
-;; extra by djseed .. psym=88 means use the already defined usersymbol
- if psym[i] eq 88 then psym[i] =8
- if vectorfont[i] ne '' then begin
-; if (num eq 1) and vertical then xp = x + xt/2 ; IF 1, CENTERED.
- xyouts,xp,yp,vectorfont[i],width=width,color=colors[i] $
- ,size=charsize,align=xalign,charthick = charthick,/norm
- xt = xt > width
- xp = xp + width/2.
- endif else begin
- if symline and (linestyle[i] ge 0) then plots,xp,yp,color=colors[i] $
- ,/normal,linestyle=linestyle[i],psym=psym[i],symsize=symsize[i], $
- thick=thick[i]
- endelse
-
- if vertical then x = x + xt else if ltor then x = max(xp) else x = min(xp)
- if symline then x = x + xspacing
- TEXT_ONLY:
- if vertical and (vectorfont[i] eq '') and symline and (linestyle[i] eq -99) then x=x0 + xspacing
- xyouts,x,y,delimiter,width=width,/norm,color=textcolors[i], $
- size=charsize,align=xalign,charthick = charthick
- x = x + width*xsign
- if width ne 0 then x = x + 0.5*xspacing
- xyouts,x,y,items[i],width=width,/norm,color=textcolors[i],size=charsize, $
- align=xalign,charthick=charthick
- x = x + width*xsign
- if not vertical and (i lt (n-1)) then x = x+2*xspacing; ADD INTER-ITEM SPACE
- xfinal = (x + xspacing*margin)
- if ltor then xend = xfinal > xend else xend = xfinal < xend ; UPDATE END X
- endfor
-
- if (iclr lt clear ) then begin
-; =====>> CLEAR AREA
- x = position[0]
- y = position[1]
- if vertical then bottom = n else bottom = 1
- ywidth = - (2*margin+bottom-0.5)*yspacing
- corners = [x,y+ywidth,xend,y]
- polyfill,[x,xend,xend,x,x],y + [0,0,ywidth,ywidth,0],/norm,color=-1
-; plots,[x,xend,xend,x,x],y + [0,0,ywidth,ywidth,0],thick=2
- endif else begin
-
-;
-; =====>> OUTPUT BORDER
-;
- x = position[0]
- y = position[1]
- if vertical then bottom = n else bottom = 1
- ywidth = - (2*margin+bottom-0.5)*yspacing
- corners = [x,y+ywidth,xend,y]
- if box then plots,[x,xend,xend,x,x],y + [0,0,ywidth,ywidth,0],/norm
- return
- endelse
-endfor
-
-end
-
diff --git a/examples/Disk-Patch/GravityOnly/run.sh b/examples/Disk-Patch/GravityOnly/run.sh
index 89e7c10ce19ce0576127551f63eea5c9e4e955ff..a123ad24d7ca34105c22f5f31e75c688c681288f 100755
--- a/examples/Disk-Patch/GravityOnly/run.sh
+++ b/examples/Disk-Patch/GravityOnly/run.sh
@@ -3,8 +3,8 @@
# 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..."
+ echo "Generating initial conditions for the disk-patch example..."
python makeIC.py 1000
fi
-../swift -g -t 2 externalGravity.yml
+../../swift -g -t 2 disk-patch.yml
diff --git a/examples/Disk-Patch/HydroStatic/readme.text b/examples/Disk-Patch/HydroStatic/README
similarity index 81%
rename from examples/Disk-Patch/HydroStatic/readme.text
rename to examples/Disk-Patch/HydroStatic/README
index 314334fd0383db0a74746da85cb3ab298aa3b060..3277115eea12885c3136d205f341086ca5d98385 100644
--- a/examples/Disk-Patch/HydroStatic/readme.text
+++ b/examples/Disk-Patch/HydroStatic/README
@@ -1,6 +1,7 @@
-genertae and evolve a disk-patch, where gas is in hydrostatic
+generate and evolve a disk-patch, where gas is in hydrostatic
equilibrium with an imposed external gravutation force, using teh
-equations from Creasey, Theuns & Bower 2013.
+equations from Creasey, Theuns & Bower, 2013, MNRAS, Volume 429,
+Issue 3, p.1922-1948
Run the swift using setting
#define EXTERNAL_POTENTIAL_DISK_PATCH
diff --git a/examples/Disk-Patch/HydroStatic/legend.pro b/examples/Disk-Patch/HydroStatic/legend.pro
deleted file mode 100755
index e510bdfcd7dcabed796cb090c6fa7d54536608b6..0000000000000000000000000000000000000000
--- a/examples/Disk-Patch/HydroStatic/legend.pro
+++ /dev/null
@@ -1,459 +0,0 @@
-;+
-; 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/legend.pro b/examples/IsothermalPotential/GravityOnly/legend.pro
deleted file mode 100755
index e510bdfcd7dcabed796cb090c6fa7d54536608b6..0000000000000000000000000000000000000000
--- a/examples/IsothermalPotential/GravityOnly/legend.pro
+++ /dev/null
@@ -1,459 +0,0 @@
-;+
-; 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/main.c b/examples/main.c
index 59f1a638284da3b3817a2fa0ab1b9e78175e3d4c..d2bf974463c2a9d44284d875160aed008f54cf07 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -329,11 +329,12 @@ int main(int argc, char *argv[]) {
/* Initialise the hydro properties */
struct hydro_props hydro_properties;
hydro_props_init(&hydro_properties, params);
- if (with_hydro && myrank == 0) eos_print();
+ if (with_hydro && myrank == 0) eos_print();
/* Initialise the external potential properties */
struct external_potential potential;
- if (with_external_gravity) potential_init(params, &prog_const, &us, &potential);
+ if (with_external_gravity)
+ potential_init(params, &prog_const, &us, &potential);
if (with_external_gravity && myrank == 0) potential_print(&potential);
/* Read particles and space information from (GADGET) ICs */
diff --git a/src/const.h b/src/const.h
index d2110be497fbe821c62911fcb3f9c3079770de0b..ec433e00a22062e887b8e720895e5f448c2b897d 100644
--- a/src/const.h
+++ b/src/const.h
@@ -48,7 +48,8 @@
/* Equation of state choice */
#define EOS_IDEAL_GAS
-/* if EOS_ISOTHERMAL_GAS is defined, keep thermal energy per unit mass equal to EOS_ISOTHERMAL_GAS in programme units */
+/* if EOS_ISOTHERMAL_GAS is defined, keep thermal energy per unit mass equal to
+ * EOS_ISOTHERMAL_GAS in programme units */
//#define EOS_ISOTHERMAL_GAS (20.2615290634)
/* Kernel function to use */
@@ -74,10 +75,10 @@
//#define EXTERNAL_POTENTIAL_POINTMASS
//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
#define EXTERNAL_POTENTIAL_DISK_PATCH
-/* Add viscuous force to gas particles to speed-up glass making for disk-patch ICs */
+/* Add viscuous force to gas particles to speed-up glass making for disk-patch
+ * ICs */
//#define EXTERNAL_POTENTIAL_DISK_PATCH_ICS
-
/* Are we debugging ? */
//#define SWIFT_DEBUG_CHECKS
diff --git a/src/equation_of_state.h b/src/equation_of_state.h
index 672b44c55d298be3408126415be90de133342656..5ca1ace6324f11b6436061dc48bc5616e393a15c 100644
--- a/src/equation_of_state.h
+++ b/src/equation_of_state.h
@@ -41,7 +41,7 @@
/* ------------------------------------------------------------------------- */
#if defined(EOS_IDEAL_GAS)
#if defined(EOS_ISOTHERMAL_GAS)
-#error: cannot have ideal and isothermal gas at the same time!
+#error : cannot have ideal and isothermal gas at the same time!
#endif
/**
@@ -128,15 +128,16 @@ gas_soundspeed_from_internal_energy(float density, float u) {
return sqrtf(u * hydro_gamma * hydro_gamma_minus_one);
}
-__attribute__((always_inline)) INLINE static void eos_print(){
- message(" assuming an equation of state of an ideal gas, with gamma = %e", hydro_gamma);
+__attribute__((always_inline)) INLINE static void eos_print() {
+ message(" assuming an equation of state of an ideal gas, with gamma = %e",
+ hydro_gamma);
}
-__attribute__((always_inline)) INLINE static float hydro_update_entropy(float density, float entropy)
-{
+__attribute__((always_inline)) INLINE static float hydro_update_entropy(
+ float density, float entropy) {
return entropy;
}
-__attribute__((always_inline)) INLINE static float hydro_update_entropy_rate(float density, float entropy)
-{
+__attribute__((always_inline)) INLINE static float hydro_update_entropy_rate(
+ float density, float entropy) {
return hydro_gamma_minus_one * pow_minus_gamma_minus_one(density);
}
/* ------------------------------------------------------------------------- */
@@ -174,7 +175,8 @@ __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
__attribute__((always_inline)) INLINE static float
gas_entropy_from_internal_energy(float density, float u) {
- return hydro_gamma_minus_one * EOS_ISOTHERMAL_GAS / pow(density, hydro_gamma-1);
+ return hydro_gamma_minus_one * EOS_ISOTHERMAL_GAS /
+ pow(density, hydro_gamma - 1);
}
/**
@@ -202,16 +204,19 @@ gas_soundspeed_from_internal_energy(float density, float u) {
error("Missing definition !");
return 0.f;
}
-__attribute__((always_inline)) INLINE static void eos_print(){
- message(" assuming an equation of state for an isothermal gas with utherm= %e and gamma = %e", EOS_ISOTHERMAL_GAS, hydro_gamma);
+__attribute__((always_inline)) INLINE static void eos_print() {
+ message(
+ " assuming an equation of state for an isothermal gas with utherm= %e "
+ "and gamma = %e",
+ EOS_ISOTHERMAL_GAS, hydro_gamma);
}
-__attribute__((always_inline)) INLINE static float hydro_update_entropy(float density, float entropy)
-{
+__attribute__((always_inline)) INLINE static float hydro_update_entropy(
+ float density, float entropy) {
float thermal_energy = EOS_ISOTHERMAL_GAS;
return gas_entropy_from_internal_energy(density, thermal_energy);
}
-__attribute__((always_inline)) INLINE static float hydro_update_entropy_rate(float density, float entropy)
-{
+__attribute__((always_inline)) INLINE static float hydro_update_entropy_rate(
+ float density, float entropy) {
return 0;
}
/* ------------------------------------------------------------------------- */
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 025c881d169186e7efd483fee0c74a601b85cc50..ca1af32096d2c392eef27cf9a9adb8c824b8ad28 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -47,7 +47,8 @@ gravity_compute_timestep_external(const struct external_potential* potential,
phys_const, gp));
#endif
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
- dt = fmin(dt, external_gravity_disk_patch_timestep(potential, phys_const, gp));
+ dt =
+ fmin(dt, external_gravity_disk_patch_timestep(potential, phys_const, gp));
#endif
return dt;
}
@@ -74,7 +75,6 @@ gravity_compute_timestep_self(const struct phys_const* const phys_const,
return FLT_MAX;
}
-
/**
* @brief Initialises the g-particles for the first time
*
@@ -117,7 +117,7 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(
* @param const_G Newton's constant
*/
__attribute__((always_inline)) INLINE static void gravity_end_force(
- struct gpart* gp, const double const_G) {
+ struct gpart* gp, const double const_G) {
/* Let's get physical... */
gp->a_grav[0] *= const_G;
@@ -135,8 +135,7 @@ __attribute__((always_inline)) INLINE static void gravity_end_force(
* @param gp The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void external_gravity(
- const double time,
- const struct external_potential* potential,
+ const double time, const struct external_potential* potential,
const struct phys_const* const phys_const, struct gpart* gp) {
#ifdef EXTERNAL_POTENTIAL_POINTMASS
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index 54f55c79aaf489942f59be2d3623627e30f672a2..b89b4f70aea1c29e2df4df3615e34552e4ea3f8d 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -21,7 +21,6 @@
#include "hydro.h"
#include "io_properties.h"
#include "kernel_hydro.h"
-#include "equation_of_state.h"
/**
* @brief Specifies which particle fields to read from a dataset
diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h
index 3145390e3fb9e3b05dadd1cac56b98cec1d59e38..0cd830c51a1836c0a56e7b4097990bf4d07a101f 100644
--- a/src/hydro/Gadget2/hydro_part.h
+++ b/src/hydro/Gadget2/hydro_part.h
@@ -64,7 +64,7 @@ struct part {
/* Derivative of the density with respect to smoothing length. */
float rho_dh;
- // union {
+ union {
struct {
@@ -100,6 +100,7 @@ struct part {
float h_dt;
} force;
+ };
/* Particle ID. */
long long id;
diff --git a/src/potentials.c b/src/potentials.c
index 26507622c3f93f86445181291172dd958f31f17a..7be7d48bc046f11179a683705a60a21e017c8a90 100644
--- a/src/potentials.c
+++ b/src/potentials.c
@@ -33,7 +33,7 @@
* @param potential The external potential properties to initialize
*/
void potential_init(const struct swift_params* parameter_file,
- const struct phys_const* const phys_const,
+ const struct phys_const* const phys_const,
struct UnitSystem* us,
struct external_potential* potential) {
@@ -67,17 +67,19 @@ void potential_init(const struct swift_params* parameter_file,
#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
- potential -> disk_patch_potential.surface_density =
- parser_get_param_double(parameter_file, "Disk-PatchPotential:surface_density");
- potential -> disk_patch_potential.scale_height =
- parser_get_param_double(parameter_file, "Disk-PatchPotential:scale_height");
- potential -> disk_patch_potential.z_disk =
- parser_get_param_double(parameter_file, "Disk-PatchPotential:z_disk");
- potential -> disk_patch_potential.timestep_mult =
- parser_get_param_double(parameter_file, "Disk-PatchPotential:timestep_mult");
- potential -> disk_patch_potential.dynamical_time = sqrt(potential->disk_patch_potential.scale_height / (phys_const->const_newton_G * potential->disk_patch_potential.surface_density));
+ potential->disk_patch_potential.surface_density = parser_get_param_double(
+ parameter_file, "Disk-PatchPotential:surface_density");
+ potential->disk_patch_potential.scale_height = parser_get_param_double(
+ parameter_file, "Disk-PatchPotential:scale_height");
+ potential->disk_patch_potential.z_disk =
+ parser_get_param_double(parameter_file, "Disk-PatchPotential:z_disk");
+ potential->disk_patch_potential.timestep_mult = parser_get_param_double(
+ parameter_file, "Disk-PatchPotential:timestep_mult");
+ potential->disk_patch_potential.dynamical_time =
+ sqrt(potential->disk_patch_potential.scale_height /
+ (phys_const->const_newton_G *
+ potential->disk_patch_potential.surface_density));
#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
-
}
/**
@@ -108,13 +110,17 @@ void potential_print(const struct external_potential* potential) {
#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
- message("Disk-patch potential properties are surface_density = %e disk height= %e scale height= %e timestep multiplier= %e",
- potential->disk_patch_potential.surface_density,
- potential->disk_patch_potential.z_disk,
- potential->disk_patch_potential.scale_height,
- potential->disk_patch_potential.timestep_mult);
+ message(
+ "Disk-patch potential properties are surface_density = %e disk height= "
+ "%e scale height= %e timestep multiplier= %e",
+ potential->disk_patch_potential.surface_density,
+ potential->disk_patch_potential.z_disk,
+ potential->disk_patch_potential.scale_height,
+ potential->disk_patch_potential.timestep_mult);
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH_ICS
- message("Disk-patch potential: imposing growth of gravity over time, and adding viscous force to gravity");
+ message(
+ "Disk-patch potential: imposing growth of gravity over time, and adding "
+ "viscous force to gravity");
#endif
#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
}
diff --git a/src/potentials.h b/src/potentials.h
index 33bc64d57f981ef8e64c84e43852fd2f77dcc4ff..422b5533c32a05065dd6c870040d96e00c82ae38 100644
--- a/src/potentials.h
+++ b/src/potentials.h
@@ -25,8 +25,8 @@
#include "../config.h"
/* Some standard headers. */
-#include <math.h>
#include <float.h>
+#include <math.h>
/* Local includes. */
#include "const.h"
@@ -35,7 +35,7 @@
#include "part.h"
#include "physical_constants.h"
#include "units.h"
-#include "math.h"
+
/* External Potential Properties */
struct external_potential {
@@ -56,11 +56,11 @@ struct external_potential {
#endif
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
struct {
- double surface_density;
- double scale_height;
- double z_disk;
- double dynamical_time;
- double timestep_mult;
+ double surface_density;
+ double scale_height;
+ double z_disk;
+ double dynamical_time;
+ double timestep_mult;
} disk_patch_potential;
#endif
};
@@ -68,17 +68,17 @@ struct external_potential {
/* external potential: disk_patch */
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
/**
- * @brief Computes the time-step due to the acceleration from a hydrostatic disk (Creasy, Theuns & Bower, 2013)
+ * @brief Computes the time-step due to the acceleration from a hydrostatic disk
+ * (Creasy, Theuns & Bower, 2013)
*
* @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(
- const struct external_potential* potential,
- const struct phys_const* const phys_const,
- const struct gpart* const g) {
+__attribute__((always_inline)) INLINE 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 = potential->disk_patch_potential.dynamical_time;
float dt = dt_dyn;
@@ -87,78 +87,85 @@ __attribute__((always_inline)) INLINE static float external_gravity_disk_patch_t
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 z_accel = 2 * M_PI * phys_const->const_newton_G *
+ potential->disk_patch_potential.surface_density *
+ 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;
- }
+ 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 dt^2 * acceleration < fraction of scale height of disk */
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);
- }
+ 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 (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.);
+ }
return potential->disk_patch_potential.timestep_mult * dt;
}
/**
- * @brief Computes the gravitational acceleration from a hydrostatic disk (Creasy, Theuns & Bower, 2013)
+ * @brief Computes the gravitational acceleration from a hydrostatic disk
+ * (Creasy, Theuns & Bower, 2013)
*
* @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 double time, const struct external_potential* potential, const struct phys_const* const phys_const, struct gpart* g) {
+__attribute__((always_inline)) INLINE static void
+external_gravity_disk_patch_potential(
+ const double time, const struct external_potential* potential,
+ const struct phys_const* const phys_const, struct gpart* g) {
const float G_newton = phys_const->const_newton_G;
- const float t_dyn = potential->disk_patch_potential.dynamical_time;
const float dz = g->x[2] - potential->disk_patch_potential.z_disk;
- g->a_grav[0] += 0;
- g->a_grav[1] += 0;
+ g->a_grav[0] += 0;
+ g->a_grav[1] += 0;
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH_ICS
float reduction_factor = 1.;
- if(time < 5 * t_dyn)
- reduction_factor = time / (5 * t_dyn);
+ if (time < 5 * t_dyn) reduction_factor = time / (5 * t_dyn);
#else
const float reduction_factor = 1.;
#endif
- const float z_accel = reduction_factor * 2 * G_newton * M_PI * potential->disk_patch_potential.surface_density
- * tanh(fabs(dz) / potential->disk_patch_potential.scale_height);
+ const float z_accel =
+ reduction_factor * 2 * G_newton * M_PI *
+ potential->disk_patch_potential.surface_density *
+ tanh(fabs(dz) / potential->disk_patch_potential.scale_height);
if (dz > 0)
- g->a_grav[2] -= z_accel / G_newton; /* returned acceleraton is multiplied by G later on */
- if (dz < 0)
- g->a_grav[2] += z_accel / G_newton;
-
+ g->a_grav[2] -=
+ z_accel /
+ G_newton; /* returned acceleraton is multiplied by G later on */
+ if (dz < 0) g->a_grav[2] += z_accel / G_newton;
- if(abs(g->id_or_neg_offset) == 1)
- message(" time= %e, rf= %e, az= %e", time, reduction_factor, g->a_grav[2]);
+ if (abs(g->id_or_neg_offset) == 1)
+ message(" time= %e, rf= %e, az= %e", time, reduction_factor, g->a_grav[2]);
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH_ICS
/* TT: add viscous drag */
- g->a_grav[0] -= g->v_full[0] / (2 * t_dyn) / G_newton;
- g->a_grav[1] -= g->v_full[1] / (2 * t_dyn) / G_newton;
- g->a_grav[2] -= g->v_full[2] / (2 * t_dyn) / G_newton;
-#endif
+ g->a_grav[0] -= g->v_full[0] / (2 * t_dyn) / G_newton;
+ g->a_grav[1] -= g->v_full[1] / (2 * t_dyn) / G_newton;
+ g->a_grav[2] -= g->v_full[2] / (2 * t_dyn) / G_newton;
+#endif
}
#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
@@ -177,129 +184,145 @@ external_gravity_isothermalpotential_timestep(
const struct phys_const* const phys_const, const struct gpart* const g) {
const struct gpart* const g) {
- const float dx = g->x[0] - potential->isothermal_potential.x;
- const float dy = g->x[1] - potential->isothermal_potential.y;
- 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];
-
- return potential->isothermal_potential.timestep_mult * sqrtf(a_2 / dota_2);
-}
-
-/**
- * @brief Computes the gravitational acceleration of a particle due to a point
- * mass
- *
- * Note that the accelerations are multiplied by Newton's G constant later on.
- *
- * @param potential The #external_potential used in the run.
- * @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_isothermalpotential(const struct external_potential* potential, const struct phys_const* const phys_const, struct gpart* g) {
-__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 G_newton = phys_const->const_newton_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;
- const double term = -vrot * vrot * rinv2 / G_newton;
-
- g->a_grav[0] += term * dx;
- g->a_grav[1] += term * dy;
- g->a_grav[2] += term * dz;
-
-}
+ 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];
+
+ return potential->isothermal_potential.timestep_mult *
+ sqrtf(a_2 / dota_2);
+ }
+
+ /**
+ * @brief Computes the gravitational acceleration of a particle due to a
+ * point
+ * mass
+ *
+ * Note that the accelerations are multiplied by Newton's G constant
+ * later on.
+ *
+ * @param potential The #external_potential used in the run.
+ * @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_isothermalpotential(
+ const struct external_potential* potential,
+ const struct phys_const* const phys_const, struct gpart* g) {
+ __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 G_newton = phys_const->const_newton_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;
+ const double term = -vrot * vrot * rinv2 / G_newton;
+
+ g->a_grav[0] += term * dx;
+ g->a_grav[1] += term * dy;
+ g->a_grav[2] += term * dz;
+ }
#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
/* Include exteral pointmass potential */
#ifdef EXTERNAL_POTENTIAL_POINTMASS
-/**
- * @brief Computes the time-step due to the acceleration from a point mass
- *
- * @param potential The properties of the externa potential.
- * @param phys_const The physical constants in internal units.
- * @param g Pointer to the g-particle data.
- */
-__attribute__((always_inline)) INLINE static float external_gravity_pointmass_timestep(const struct external_potential* potential,
- const struct phys_const* const phys_const,
- const struct gpart* const g) {
-
- const float G_newton = phys_const->const_newton_G;
- const float dx = g->x[0] - potential->point_mass.x;
- const float dy = g->x[1] - potential->point_mass.y;
- const float dz = g->x[2] - potential->point_mass.z;
- const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
- const float drdv = (g->x[0] - potential->point_mass.x) * (g->v_full[0]) +
- (g->x[1] - potential->point_mass.y) * (g->v_full[1]) +
- (g->x[2] - potential->point_mass.z) * (g->v_full[2]);
- const float dota_x = G_newton * potential->point_mass.mass * rinv * rinv *
- rinv * (-g->v_full[0] + 3.f * rinv * rinv * drdv * dx);
- const float dota_y = G_newton * potential->point_mass.mass * rinv * rinv *
- rinv * (-g->v_full[1] + 3.f * rinv * rinv * drdv * dy);
- const float dota_z = G_newton * potential->point_mass.mass * rinv * rinv *
- rinv * (-g->v_full[2] + 3.f * rinv * rinv * drdv * dz);
- 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];
-
- return potential->point_mass.timestep_mult * sqrtf(a_2 / dota_2);
-}
-
-/**
- * @brief Computes the gravitational acceleration of a particle due to a point
- * mass
- *
- * Note that the accelerations are multiplied by Newton's G constant later on.
- *
- * @param potential The proerties of the external 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(
- const struct external_potential* potential,
- const struct phys_const* const phys_const, struct gpart* g) {
-
- const float dx = g->x[0] - potential->point_mass.x;
- const float dy = g->x[1] - potential->point_mass.y;
- const float dz = g->x[2] - potential->point_mass.z;
- const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
- const float rinv3 = rinv * rinv * rinv;
+ /**
+ * @brief Computes the time-step due to the acceleration from a point mass
+ *
+ * @param potential The properties of the externa potential.
+ * @param phys_const The physical constants in internal units.
+ * @param g Pointer to the g-particle data.
+ */
+ __attribute__((always_inline)) INLINE static float
+ external_gravity_pointmass_timestep(
+ const struct external_potential* potential,
+ const struct phys_const* const phys_const,
+ const struct gpart* const g) {
- g->a_grav[0] += -potential->point_mass.mass * dx * rinv3;
- g->a_grav[1] += -potential->point_mass.mass * dy * rinv3;
- g->a_grav[2] += -potential->point_mass.mass * dz * rinv3;
-}
+ const float G_newton = phys_const->const_newton_G;
+ const float dx = g->x[0] - potential->point_mass.x;
+ const float dy = g->x[1] - potential->point_mass.y;
+ const float dz = g->x[2] - potential->point_mass.z;
+ const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
+ const float drdv = (g->x[0] - potential->point_mass.x) * (g->v_full[0]) +
+ (g->x[1] - potential->point_mass.y) * (g->v_full[1]) +
+ (g->x[2] - potential->point_mass.z) * (g->v_full[2]);
+ const float dota_x = G_newton * potential->point_mass.mass * rinv * rinv *
+ rinv *
+ (-g->v_full[0] + 3.f * rinv * rinv * drdv * dx);
+ const float dota_y = G_newton * potential->point_mass.mass * rinv * rinv *
+ rinv *
+ (-g->v_full[1] + 3.f * rinv * rinv * drdv * dy);
+ const float dota_z = G_newton * potential->point_mass.mass * rinv * rinv *
+ rinv *
+ (-g->v_full[2] + 3.f * rinv * rinv * drdv * dz);
+ 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];
+
+ return potential->point_mass.timestep_mult * sqrtf(a_2 / dota_2);
+ }
+
+ /**
+ * @brief Computes the gravitational acceleration of a particle due to a
+ * point
+ * mass
+ *
+ * Note that the accelerations are multiplied by Newton's G constant later
+ * on.
+ *
+ * @param potential The proerties of the external 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(const struct external_potential* potential,
+ const struct phys_const* const phys_const,
+ struct gpart* g) {
+
+ const float dx = g->x[0] - potential->point_mass.x;
+ const float dy = g->x[1] - potential->point_mass.y;
+ const float dz = g->x[2] - potential->point_mass.z;
+ const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
+ const float rinv3 = rinv * rinv * rinv;
+
+ g->a_grav[0] += -potential->point_mass.mass * dx * rinv3;
+ g->a_grav[1] += -potential->point_mass.mass * dy * rinv3;
+ g->a_grav[2] += -potential->point_mass.mass * dz * rinv3;
+ }
#endif /* EXTERNAL_POTENTIAL_POINTMASS */
-/* Now, some generic functions, defined in the source file */
-void potential_init(const struct swift_params* parameter_file,
- const struct phys_const* const phys_const,
- struct UnitSystem* us,
- struct external_potential* potential);
+ /* Now, some generic functions, defined in the source file */
+ void potential_init(const struct swift_params* parameter_file,
+ const struct phys_const* const phys_const,
+ struct UnitSystem* us,
+ struct external_potential* potential);
-void potential_print(const struct external_potential* potential);
+ void potential_print(const struct external_potential* potential);
#endif /* SWIFT_POTENTIALS_H */
diff --git a/src/timestep.h b/src/timestep.h
index fc1869156fc356f51b924379b5f4666233ab0491..cb72134f54cfc3e3465796311d3b6328472dd8c6 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -112,12 +112,11 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
new_dt_grav = fminf(new_dt_external, new_dt_self);
- if(p->id == -1)
- {
- printParticle_single(p, xp);
- message(" dt_hydro= %e, dt_grav_external= %e dt_grav_self=%e ",new_dt_hydro,new_dt_external,new_dt_self);
- }
-
+ if (p->id == -1) {
+ printParticle_single(p, xp);
+ message(" dt_hydro= %e, dt_grav_external= %e dt_grav_self=%e ",
+ new_dt_hydro, new_dt_external, new_dt_self);
+ }
}
/* Final time-step is minimum of hydro and gravity */
@@ -130,9 +129,7 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
: FLT_MAX;
new_dt = fminf(new_dt, dt_h_change);
- if(p->id == -1)
- message(" new_dt= %e", new_dt);
-
+ if (p->id == -1) message(" new_dt= %e", new_dt);
/* Limit timestep within the allowed range */
new_dt = fminf(new_dt, e->dt_max);