comp.lang.idl-pvwave archive: archive » Gradient of two dimensional field

Home » Public Forums » archive » Gradient of two dimensional field

Show: Today's Messages :: Show Polls :: Message Navigator
E-mail to friend

Re: Gradient of two dimensional field [message #8293 is a reply to message #8200]

Wed, 19 February 1997 00:00

Andy Loughe
Messages: 174
Registered: November 1995

Senior Member

David Fanning wrote:
>
> Andy Loughe <afl@cdc.noaa.gov> writes in response to Wilpert Martin:
>
>>> we want to determine the electrical field from a given potential,
>>> i.e. we have to calculate the gradient of a two dimensional array.
>
>> I would think that the shift function (used twice)
>> could be used to do this.
>
> Andy, do you think you could you give those of us who are
> wondering about this just a small example of what you
> mean? Thanks!
>
> David

A simple solution was posted already, but I think it accomplished the
task via one-sided differences. To use two-sided differences one may
wish to use the shift function twice for each vector component, then
use shift again to handle the pesky boundaries (assuming non-cyclic
boundary values), then compute the "magnitude" of the gradient, and
return this value.

Here is a soulution that I have not had time to debug.
If improvements need to be made, please let me know...

; Compute the vector magnitude of the gradient.
;
; Andrew F. Loughe (afl@cdc.noaa.gov)
;

function grad, data, x, y

sz = size(data) & im = sz(1) & jm = sz(2)

if ( N_params() eq 0 ) then message, ' grad_data = grad(data, x, y)'
if ( sz(0) ne 2 ) then message, 'Input data must be 2-D'
if ( N_elements(x) eq 0) then x = indgen(im) + 1
if ( N_elements(y) eq 0) then y = jm - indgen(jm)

; Begin here

x2 = x # replicate(1, jm)
y2 = replicate(1,jm) # y

; i-component
if (x(1) gt x(0)) then dx = shift(x2, -1, 0) - shift(x2, 1, 0)
if (x(1) lt x(0)) then dx = shift(x2, 1, 0) - shift(x2, -1, 0)
grad_x = ( shift(data, -1, 0) - shift(data, 1, 0) ) / dx

; j-component
if (y(1) lt y(0)) then dy = shift(y2, 0, 1) - shift(y2, 0, -1)
if (y(1) gt y(0)) then dy = shift(y2, 0, -1) - shift(y2, 0, 1)
grad_y = ( shift(data, 0, 1) - shift(data, 0, -1) ) / dy

; But for non-cyclic boundary values we still have a problem...
; Take care of the outer rows and columns
grad_y(*,jm-1) = ( data(*,jm-1) - data(*,jm-2) ) / $
( y2(*,jm-1) - y2(*,jm-2) )
grad_y(*,0) = ( data(*,1) - data(*,0) ) / ( y2(*,1) - y2(*,0) )

grad_x(im-1,*) = ( data(im-1,*) - data(im-2,*) ) / $
( x2(im-1,*) - x2(im-2,*) )
grad_x(0,*) = ( data(1,*) - data(0,*) ) / ( x2(1,*) - x2(0,*) )

grad = sqrt(grad_x^2 + grad_y^2)

return, grad

end

Report message to a moderator

[Message index]

		Gradient of two dimensional field By: Wilpert_Martin on Wed, 19 February 1997 00:00
		Re: Gradient of two dimensional field By: Andy Loughe on Wed, 19 February 1997 00:00
		Re: Gradient of two dimensional field By: davidf on Wed, 19 February 1997 00:00
		Re: Gradient of two dimensional field By: stockwell on Wed, 19 February 1997 00:00
		Re: Gradient of two dimensional field By: Andy Loughe on Wed, 19 February 1997 00:00
		Re: Gradient of two dimensional field By: bowman on Wed, 19 February 1997 00:00
		Re: Gradient of two dimensional field By: brian.jackel on Wed, 19 February 1997 00:00
		Re: Gradient of two dimensional field By: Mirko Vukovic on Wed, 19 February 1997 00:00
		Re: Gradient of two dimensional field By: Andy Loughe on Mon, 03 March 1997 00:00

Previous Topic:	Re: Printing array on one line
Next Topic:	using GET_KBRD

-=] Back to Top [=-

[ Syndicate this forum (XML) ] [

] [

]

Current Time: Tue Dec 02 11:11:56 PST 2025

Total time taken to generate the page: 1.04571 seconds