Re: Accuracy problem [message #47948 is a reply to message #19921] |
Fri, 10 March 2006 22:52   |
mmeron
Messages: 44 Registered: October 2003
|
Member |
|
|
In article <MPG.1e7c116c2aa1c294989bd5@news.frii.com>, David Fanning <davidf@dfanning.com> writes:
> Juan Arrieta writes:
>
>> I am preparing a simple code to transform between cartesian vectors
>> and Keplerian elements (semimajor axis, inclination, eccentricity, and
>> so forth). This is a simple problem in astrodynamics.
>>
>> At some point in my code, I need to obtain unit vectors. For instance,
>> a line of the code is:
>>
>> U0 = ACOS( (TRANSPOSE(NDVCT) # R) / ( NORM(NDVCT) * NORM(R) ) )
>> % Program caused arithmetic error: Floating illegal operand
>> print,U0
>> NaN
>>
>> Any comments as for what would the problem be? This seems like a
>> roundoff error somewhere in the program, but I am not doing anything
>> "fancy" here.
>
> Getting errors with computers (alas!) doesn't require much
> "fancy" programming. I don't know exactly what the problem
> is here, but clearly you need more precision. I'd start by
> doing everything in DOUBLE precision. You don't tell us
> what datatype NDVCT and R are but the NORM function is
> being done as a FLOAT, since you haven't set the DOUBLE
> keyword.
>
> I'd step back a couple of steps before you introduced us
> to the problem and make sure everything is done in double
> precision values. Then, after you have convinced yourself
> you've done everything humanly possible, I think you might
> be justified in confining the arguments to ACOS to the
> proper range:
>
> arg = ACOS(0.0 > expr < 1.0)
>
> Cheers,
>
> David
>
> P.S. I will say that I run into a LOT of floating underflow
> errors when I am working with astronomy data. Mostly when
> values get close to 0. (I see a -0.000000 value in your
> example.) Maybe these are generated when images are flat-fielded,
> or processed in other ways. I've even resorted to cleaning
> these "close to zero" values up before I work with the images.
> It tends to keep the blood pressure down a little.
>
> I = Where(image GT -1e-8 AND image LT 1e-8, count)
> if count GT 0 THEN image[I] = 0.0
>
> I just checked my pulse, and look at that, my blood pressure
> is rising just *thinking* of those damn floating underflow messages!
> --
Well, they're a nuisance, or rather used to be one. I've a little
routine in my library (called FPU_FIX) which is killing them on sight
so nowadays I hardly ever see them.
Mati Meron | "When you argue with a fool,
meron@cars.uchicago.edu | chances are he is doing just the same"
|
|
|