;+
; NAME:
;       POTENTIAL_NEWTON
;
; PURPOSE:
;       Calculates the gravitational potential from a particle distribution
;       at a list of positions.
;
; CATEGORY:
;       Astro
;
; CALLING SEQUENCE:
;       Result = POTENTIAL_NEWTON(X, Y, Z, Xpart, Ypart, Zpart, Mass)
;
; INPUTS:
;       X:     X coordinates at which to calculate potential.
;
;       Y:     Y coordinates at which to calculate potential.
;
;       Z:     Z coordinates at which to calculate potential.
;
;       Xpart: X coordinates of particle positions defining the mass
;              distribution.
;
;       Ypart: Y coordinates of particle positions defining the mass
;              distribution.
;
;       Zpart: Z coordinates of particle positions defining the mass
;              distribution.
;
;       Mass:  Mass of each particle.
;
; KEYWORD PARAMETERS:
;       LENGTHUNIT: Length unit, in cm (or kpc if /ASTRO is set). Default: 1kpc.
;
;       MASSUNIT:   Mass unit, in grams (or solar masses if /ASTRO is
;                   set). Default: 1 solar mass.
;
;       SOFTENING:  Plummer softening length. Default: 0 (no softening).
;
;       ASTRO:      If /ASTRO is set then LENGTHUNIT and MASSUNIT are given
;                   in kpc and solar masses respectively, otherwise they are
;                   in CGS (Lengthunit in cm, Massunit in grams).
;
;       LOMEM:      If /LOMEM is set then sacrifice efficiency for
;                   lower memory usage
;
; OUTPUTS:
;       This function returns the gravitational potential at each X,Y,Z position.
;       If there are NPOS positions, Result is a vector of length NPOS.
;
; EXAMPLE:
;       Calculate the gravitational potential on a line along the x-axis from
;       6 point masses placed at the vertices of a cube:
;
;       xmasspos = [-1,-1,-1,-1,1,1,1,1]
;       ymasspos = [-1,-1,1,1,-1,-1,1,1]
;       zmasspos = [-1,1,-1,1,-1,1,-1,1]
;       masses = REPLICATE(1.,8)
;       xlinepos = 0.1*FINDGEN(20)
;       ylinepos = REPLICATE(0.,20)
;       zlinepos = REPLICATE(0.,20)
;       potentials = POTENTIAL_NEWTON(xlinepos, ylinepos, zlinepos, xmasspos, ymasspos,
;         zmasspos, masses)
;
; MODIFICATION HISTORY:
;       Written by:    Jeremy Bailin
;       10 June 2008   Public release in JBIU
;
;-
function potential_newton, x, y, z, xpart, ypart, zpart, mass, $
  lengthunit=lengthunit, massunit=massunit, astro=astro_units, softening=eps, $
  lomem=lomemp

; copied from astroconst.idl to reduce dependencies:
A_pc = 3.08567758d18
A_msun = 1.9889d33
A_G = 6.673d-8

npart = n_elements(xpart)
npos = n_elements(x)

if n_elements(y) ne npos then message, 'X, Y, and Z must have the same number of elements.'
if n_elements(z) ne npos then message, 'X, Y, and Z must have the same number of elements.'
if n_elements(ypart) ne npart then message, 'XPART, YPART, and ZPART must have the same number of elements.'
if n_elements(zpart) ne npart then message, 'XPART, YPART, and ZPART must have the same number of elements.'
if n_elements(mass) ne npart then message, 'MASS must have the same number of elements as XPART.'



if n_elements(lengthunit) eq 0 then lengthunit=1e3*A_pc else begin
  if keyword_set(astro) then lengthunit*=1e3*A_pc
endelse
if n_elements(massunit) eq 0 then massunit=A_msun else begin
  if keyword_set(astro) then massunit*=A_msun
endelse
lengthunit=double(lengthunit)
massunit=double(massunit)


if n_elements(eps) eq 0 then eps=0
eps *= lengthunit


if keyword_set(lomemp) then begin
  potential = dblarr(npos)
  ; loop over whichever has fewer elements: positions or particles
  if npos lt npart then begin
    for i=0L,npos-1 do begin
      xdiff = double(x[i]*lengthunit - xpart*lengthunit)
      ydiff = double(y[i]*lengthunit - ypart*lengthunit)
      zdiff = double(z[i]*lengthunit - zpart*lengthunit)
      rdiff = sqrt(xdiff^2 + ydiff^2 + zdiff^2)+eps

      potential[i] = -A_G * total(mass*massunit / rdiff)
   endfor
  endif else begin
    for i=0L,npart-1 do begin
      xdiff = double(x*lengthunit - xpart[i]*lengthunit)
      ydiff = double(y*lengthunit - ypart[i]*lengthunit)
      zdiff = double(z*lengthunit - zpart[i]*lengthunit)
      rdiff = sqrt(xdiff^2 + ydiff^2 + zdiff^2)+eps

      potential += -A_G * mass[i]*massunit / rdiff
    endfor
  endelse

endif else begin  ; use as much memory as you want!
  xdiff = double(rebin(reform(x*lengthunit,npos,1),npos,npart) - $
    rebin(reform(xpart*lengthunit,1,npart),npos,npart))
  ydiff = double(rebin(reform(y*lengthunit,npos,1),npos,npart) - $
    rebin(reform(ypart*lengthunit,1,npart),npos,npart))
  zdiff = double(rebin(reform(z*lengthunit,npos,1),npos,npart) - $
    rebin(reform(zpart*lengthunit,1,npart),npos,npart))
  rdiff = sqrt(xdiff^2 + ydiff^2 + zdiff^2)+eps
  
  potential = -A_G * total(rebin(reform(mass*massunit,1,npart),npos,npart) / $
    rdiff, 2)
endelse

return, potential

end