Commit 96d724b8 authored by thomas.forbriger's avatar thomas.forbriger

[CLOSE] (issue7): issue7 is solved

see comment to commit 33d32e44
parent 383e7c00
this is <issue7_refract_catch_zero_offset.txt>
============================================================================
refract: catch zero offset with power-law scaling
-------------------------------------------------
Refract offers a power-law amplitude scaling. This must fail in cases where
one of the traces is located at zero offset. In such cases the entire plot
fails, which can confuse users not aware of this side-effect. Catching such
situations by a proper warning or error message might be helpful.
Options to handle the problem:
- output warning and abort
- output warning only
- forced adjustment of scaling factor for near-offset traces to an upper limit
given as a multiple of the smallest scaling factor
- elimination of near-offset traces from set of traces to be plot
- refer to second nearest offset trace
Scaling involves the following functions:
- refract_mpc.f provides subroutine mpcfactors
- refract_settracevp.f provides function settracevp
Scaling is done by adjusting the viewport for each trace individually, not by
multiplication of the samples with the mpc factors. For this reason for any
exponent (except exponent being exactly zero) a division by zero takes place
for zero-offset traces.
The exception might by just by chance. How is zero to the power of zero
defined?
From the man page of the pow-function:
- The pow() function returns the value of x raised to the power of y.
- If y is 0, the result is 1.0 (even if x is a NaN).
We potentially face three problems:
1 When using the zero offset (least offset) trace as a reference for offset
dependent scaling (as is the default), a division by zero happens in
subroutine mpcfactors.
2 Zero offset traces may come with an mpc-scaling factor equal to zero. This
means, the trace should be down-scaled to a straight line in any case.
The current concept of trace scaling does this by leaving the sample values
of the trace as they are, by using a viewport for the trace according to
current clipping (derived in function settracevp) and by setting the world
coordinates for the trace's viewport appropriately (i.e. to infinity). In
the latter step there happens the nasty division by zero (division by mpc
factor). This problem is present, even is the zero offset trace is not used
as a reference.
3 For a negative exponent the mpc factor of the zero offset trace would become
infinite even for a finite reference value.
Solutions:
1 The first problem is resolved by not using a zero-offset trace as a
reference. I recommend to make the method of taking the average amplitude in
an offset range the standard in mode 2 and 3. The default range could be
offset zero to a finite upper offset value which must not be smaller then 5
per cent of the largest offset.
2 A straight line should be plotted actually. The most straight forward
solution would be to set world coordinate range for the trace window to the
largest possible value in settracevp.
3 The zero offset trace should be removed from the set of plottable traces.
This can be done by simply returning .false. from settracevp
Action
------
- Problem 1 solved by referring amplitude scaling factors to the average
scaled amplitude in a given offset range in scaling mode 2 and 3
- Problem 2 is solved by setting an upper limit to the world coordinate range
(in counts) in the trace's viewport, such that the trace is displayed as a
straight line without causing numerical overflow.
- Problem 3 solved by enforcing non-negative exponent value
Note
----
Fortran 95 and later offers a function huge().
See gfortran documentation.
The return value can be used to adjust world coordinate limits.
`HUGE(X)' returns the largest number that is not an infinity in
the model of the type of `X'.
huge(0.)=3.40282347E+38
huge(0.d0)=1.7976931348623157E+308
settracevp make calculations of the following form:
a=b/c
where c might be too small such that a becomes NaN.
Replace this by
if abs(huge(0.)*c) < abs(b)
a=sign(huge(0.),b)
else
a=b/c
endif
a is calculated for upper an lower bound.
extend caluculation should ensure a_upper/a_lower=b_upper/b_lower
----- END OF issue7_refract_catch_zero_offset.txt -----
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment