;+ ; NAME: ; FORCE_NEWTON ; ; PURPOSE: ; Calculates the gravitational force from a particle distribution ; at a list of positions. ; ; CATEGORY: ; Astro ; ; CALLING SEQUENCE: ; Result = FORCE_NEWTON(X, Y, Z, Xpart, Ypart, Zpart, Mass) ; ; INPUTS: ; X: X coordinates at which to calculate forces. ; ; Y: Y coordinates at which to calculate forces. ; ; Z: Z coordinates at which to calculate forces. ; ; 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). ; ; OUTPUTS: ; This function returns the gravitational force at each X,Y,Z position. ; If there are NPOS positions, Result is an NPOSx3 matrix. ; ; EXAMPLE: ; Calculate the gravitational force 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) ; forces = FORCE_NEWTON(xlinepos, ylinepos, zlinepos, xmasspos, ymasspos, ; zmasspos, masses) ; ; MODIFICATION HISTORY: ; Written by: Jeremy Bailin ; 10 June 2008 Public release in JBIU ; ;- function force_newton, x, y, z, xpart, ypart, zpart, mass, $ lengthunit=lengthunit, massunit=massunit, astro=astro_units, softening=given_eps ; 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 if n_elements(given_eps) eq 0 then eps=0. else eps=double(given_eps) eps *= lengthunit 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) notitself = where(rdiff ne 0.) xdiff[notitself] /= rdiff[notitself] ydiff[notitself] /= rdiff[notitself] zdiff[notitself] /= rdiff[notitself] rdiff2 = reform(temporary(rdiff^2) + eps^2,npos,npart) forcex = - A_G * total(reform(xdiff,npos,npart) * reform(rebin(reform(mass*massunit,1,npart),npos,npart),npos,npart) / rdiff2, 2) forcey = - A_G * total(reform(ydiff,npos,npart) * reform(rebin(reform(mass*massunit,1,npart),npos,npart),npos,npart) / rdiff2, 2) forcez = - A_G * total(reform(zdiff,npos,npart) * reform(rebin(reform(mass*massunit,1,npart),npos,npart),npos,npart) / rdiff2, 2) return, [[forcex],[forcey],[forcez]] end