| Re: Force Field [message #45491] |
Wed, 14 September 2005 04:03 |
panblosky
Messages: 17
Registered: May 2005
|
Junior Member |
|
|
>
> Assuming that you are happy with first-order, centered differences, you
> want something like this
>
> grad_x = (SHIFT(phi, -1, 0) - SHIFT(phi, 1, 0))/(2.0*dx)
> grad_y = (SHIFT(phi, 0, -1) - SHIFT(phi, 0, 1))/(2.0*dy)
>
> where phi is the potential It think this handles the periodic boundary
> conditions properly, but check to make sure.
>
> Cheers, Ken Bowman
Thanks Ken, it does work, and really fast!!
|
|
|
|
| Re: Force Field [message #45494 is a reply to message #45491] |
Tue, 13 September 2005 19:13  |
Kenneth P. Bowman
Messages: 585
Registered: May 2000
|
Senior Member |
|
|
In article <1126634205.227520.96190@f14g2000cwb.googlegroups.com>,
"Andres" <panblosky@gmail.com> wrote:
> Hi all,
>
> I have a potential field on a grid (it has periodic boundary
> conditions), and I want to calculate the acceleration (force field) of
> this, i.e., f=-grad(potential), I need the gradient of the potential .
> I wrote a procedure, but it has to many for loops in it. Does any body
> know a fast and smart way to do this?
> Thanks!!
>
> Andres
Assuming that you are happy with first-order, centered differences, you
want something like this
grad_x = (SHIFT(phi, -1, 0) - SHIFT(phi, 1, 0))/(2.0*dx)
grad_y = (SHIFT(phi, 0, -1) - SHIFT(phi, 0, 1))/(2.0*dy)
where phi is the potential It think this handles the periodic boundary
conditions properly, but check to make sure.
Cheers, Ken Bowman
|
|
|
|
| Re: Force Field [message #45498 is a reply to message #45494] |
Tue, 13 September 2005 11:22 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
Andres writes:
> I have a potential field on a grid (it has periodic boundary
> conditions), and I want to calculate the acceleration (force field) of
> this, i.e., f=-grad(potential), I need the gradient of the potential .
> I wrote a procedure, but it has to many for loops in it. Does any body
> know a fast and smart way to do this?
Here is a function I use to calculate the gradient of a 2D image:
FUNCTION Gradient, image
; gradient ~= |Gx| + |Gy|
COMPILE_OPT idl2
; Check parameters.
ON_ERROR, 2
IF N_Elements(image) EQ 0 THEN Message, 'Must pass 2D image data.'
; Calculate gradient Gx.
xkernel = [ [-1.0, 0.0, 1.0], [-2.0, 0.0, 2.0], [-1.0, 0.0, 1.0] ]
gx = Convol( Float(image), xkernel, Center=1, /Edge_Wrap )
; Calculate gradient Gy.
ykernel = [ [-1.0, -2.0, -1.0], [ 0.0, 0.0, 0.0], [ 1.0, 2.0, 1.0] ]
gy = Convol( Float(image), ykernel, Center=1, /Edge_Wrap )
RETURN, ABS(gx) + ABS(gy)
END
I use it like this:
IDL> TVSCL, Gradient(image2d)
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
