| Looking for more ideas on code ... [message #32276] |
Mon, 30 September 2002 13:26  |
jeyadev
Messages: 78
Registered: February 1995
|
Member |
|
|
I have a question about how best (style and function, if possible!) to
write code for a function that has limits that have to be treated in a
special way. Consider the function
f(x) = sin(x)/x
as an example. Now, if x is always a scalar, then on just tests to see
if it is zero, and then handle that special case using a if .. then ..
else construct. But, what if x can also be scalar? I have the following
code that works:
------------------------------------------
function sinc, y
if(n_elements(y) eq 1) then begin ; y is a scalar
if(y eq 0.0) then profile = 1.0 else begin
profile = sin(y)/y
endelse
endif else begin ; y is a vector
zeros = where(y eq 0.0, ind)
if(ind gt 0) then y(zeros) = 1.0e-10 ; set zeroes to a small quantity
profile = sin(y)/y
endelse
profile = profile*profile/a0
return, profile
end
------------------------------------------
I guess the one can always set
profile(zeros) = 1.0
to handle the more general cases. But, the real question is there
a better way than
zeros = where(y eq 0.0, ind)
if(ind gt 0) then y(zeros) = "special values"
notzeros = where(y ne 0.0, ind)
if(ind gt 0) then y(notzeros) = "general definition"
I do understand that one should not compare reals, etc., but I will
clean up the numerics later.
thanks
sj
--
Surendar Jeyadev jeyadev@wrc.xerox.bounceback.com
Remove 'bounceback' for email address
|
|
|
|
|
|
|
|
| Re: Looking for more ideas on code ... [message #32338 is a reply to message #32276] |
Wed, 02 October 2002 08:26  |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
Craig Markwardt (craigmnet@cow.physics.wisc.edu) writes:
> Right, and I *always* check the count returned from WHERE, except, it
> seems, when posting to the newsgroup... Thanks for catching that.
Uh, I always thought the rule was that anything we post
to the newsgroup is *pseudo* code. :-(
Cheers,
David
--
David W. Fanning, Ph.D.
Fanning Software Consulting, Inc.
Phone: 970-221-0438, E-mail: david@dfanning.com
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Toll-Free IDL Book Orders: 1-888-461-0155
|
|
|
|
| Re: Looking for more ideas on code ... [message #32340 is a reply to message #32276] |
Wed, 02 October 2002 08:13 |
Craig Markwardt
Messages: 1869
Registered: November 1996
|
Senior Member |
|
|
thompson@orpheus.nascom.nasa.gov (William Thompson) writes:
[ Markwardt: ]
>>> profile(wh) = sin(y(wh))/y(wh)
>
[ Jeyadev : ]
>> Nice one that, when the only exceptional value is the same for
>> all "problem" points.
>
> The only suggestion I would make would be to change the last line to
>
> if ct gt 0 then profile(wh) = sin(y(wh))/y(wh)
>
> I've been bitten by that one on countless occasions.
Right, and I *always* check the count returned from WHERE, except, it
seems, when posting to the newsgroup... Thanks for catching that.
Craig
--
------------------------------------------------------------ --------------
Craig B. Markwardt, Ph.D. EMAIL: craigmnet@cow.physics.wisc.edu
Astrophysics, IDL, Finance, Derivatives | Remove "net" for better response
------------------------------------------------------------ --------------
|
|
|
|
| Re: Looking for more ideas on code ... [message #32345 is a reply to message #32276] |
Wed, 02 October 2002 07:24 |
thompson
Messages: 584
Registered: August 1991
|
Senior Member |
|
|
jeyadev@wrc.xerox.bounceback.com (Surendar Jeyadev) writes:
> In article <onr8fbugea.fsf@cow.physics.wisc.edu>,
> Craig Markwardt <craigmnet@cow.physics.wisc.edu> wrote:
>>
>> jeyadev@wrc.xerox.bounceback.com (Surendar Jeyadev) writes:
>>
>>> ....
>>>
>>> ------------------------------------------
>>>
>>> function sinc, y
>>>
>>>
>>> if(n_elements(y) eq 1) then begin ; y is a scalar
>>> if(y eq 0.0) then profile = 1.0 else begin
>>> profile = sin(y)/y
>>> endelse
>>> endif else begin ; y is a vector
>>> zeros = where(y eq 0.0, ind)
>>> if(ind gt 0) then y(zeros) = 1.0e-10 ; set zeroes to a small quantity
>>> profile = sin(y)/y
>>> endelse
>>>
>>> profile = profile*profile/a0
> My mistake there. I am actually after sinc^2 ..... but it
> has not caused any harm!
>
>>> return, profile
>>>
>>> end
>>>
>>
>> Second of all, you can simplify your logic a little, by pre-filling
>> the array with the "special case:"
>>
>> profile = y*0 + 1. ;; Tricky way to get array filled with zeroes
>> wh = where(y NE 0, ct)
>> profile(wh) = sin(y(wh))/y(wh)
> Nice one that, when the only exceptional value is the same for
> all "problem" points.
The only suggestion I would make would be to change the last line to
if ct gt 0 then profile(wh) = sin(y(wh))/y(wh)
I've been bitten by that one on countless occasions.
William Thompson
|
|
|
|
| Re: Looking for more ideas on code ... [message #32346 is a reply to message #32276] |
Wed, 02 October 2002 07:13 |
meinel
Messages: 14
Registered: February 1994
|
Junior Member |
|
|
jeyadev@wrc.xerox.bounceback.com (Surendar Jeyadev) wrote in message news:<anac2a$53h$1@news.wrc.xerox.com>...
> function sinc, y
>
>
> if(n_elements(y) eq 1) then begin ; y is a scalar
> if(y eq 0.0) then profile = 1.0 else begin
> profile = sin(y)/y
> endelse
> endif else begin ; y is a vector
> zeros = where(y eq 0.0, ind)
> if(ind gt 0) then y(zeros) = 1.0e-10 ; set zeroes to a small quantity
> profile = sin(y)/y
> endelse
>
> profile = profile*profile/a0
>
>
> return, profile
>
> end
Neat tricks! But, how about the one-line solution:
profile = (sin(y) + (y EQ 0))/(y + (y EQ 0))
Ed M
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
