Messages from Usenet group comp.lang.idl-pvwave, compiled by Paulo Penteado
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/
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
> > 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!!