Commit 33d32e44 authored by thomas.forbriger's avatar thomas.forbriger

[MERGE] (master|issue7_refract): solve zero-offset scaling issues in refract

Merge branch 'issue7_refract'

Advance refract version number to 4.14

This merge commit solves issues when scaling traces with an offset-dependent
scaling factor. The modifications are prepared on branch issue7_refract.
The issue observed with the previous version of refract was:

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.

This problem actually was caused by three problems in the code:
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 if 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.

These problems are addressed in the following way be this merge commit:
1 Problem is solved by referring amplitude scaling factors to the average
  scaled amplitude in a given offset range in scaling mode 2 and 3 (as was
  provided only in mode 3 so far), which now is the default way, the scaling
  reference is defined.
2 Problem 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. This is achieved by
  applying a lower limit ti mpc factors in subroutine refract_settracevp
3 Problem is solved by enforcing non-negative exponent value.
parents d11ca72c 5f44aa20
......@@ -90,7 +90,8 @@ REFSUB=refract_readdata.o refract_skipdata.o refract_setdefaults.o \
refract_varplot.o refract_pgparameters.o refract_ttreduce.o \
refract_message.o refract_dopicks.o refract_domodel.o refract_usage.o \
refract_cmdopt.o refract_selfilestyle.o refract_preread.o \
refract_pgdefaults.o refract_vpframe.o refract_plotoffset.o
refract_pgdefaults.o refract_vpframe.o refract_plotoffset.o \
refract_setrefrange.o refract_avgamp.o
REFOBS=refract.o $(addprefix sub/, $(REFSUB))
# Fortran dependencies
......
......@@ -30,6 +30,8 @@ c 14/11/2011 V1.4 remember whether minoff is forced on command line
c 12/11/2012 V1.5 optionally do not time shift offset shifted traces
c 12/11/2012 V1.6 added frame parameter
c 24/10/2013 V1.7 added parameters for scaling mode -S3
c 07/12/2015 V1.8 refering amplitude scaling to average in offest range
c becomes the standard (without alternative)
c
c==============================================================================
c
......@@ -54,8 +56,8 @@ c minimum offset step to be taken as two different positions
c true if plpar_minoff is forced by command line
logical plpar_forceminoff
c
logical plflag_m3avg
real plpar_m3avgxmin, plpar_m3avgxmax
c offset range to calculate reference amplitude for amplitude scaling
real plpar_avgrefxmin, plpar_avgrefxmax
c
c control modes
c -------------
......@@ -160,7 +162,7 @@ c common blocks
& plflag_seistyle, plflag_ttstyle, plpar_radius,
& plflag_hypoffset, plflag_tracenum,
& plflag_tracename, plpar_forceminoff,
& plflag_m3avg, plpar_m3avgxmin, plpar_m3avgxmax
& plpar_avgrefxmin, plpar_avgrefxmax
common /refract_elem/ elem_modbox, elem_filenames, elem_version,
& elem_annot, elem_scales, elem_data,
& elem_syntt, elem_picks, plflag_subscale,
......
......@@ -52,6 +52,13 @@ c unusual cases (see refract_setscale.f)
c 20/11/2012 V4.9 several new plot style options are implemented
c 24/10/2013 V4.10 added alternative definitions of ordinate scale
c 21/03/2014 V4.12 optionally reverse order of legend strings (-TR)
c 07/12/2015 V4.13 add option -SN; make amplitude scaling with
c respect to offset range the standard (without
c alternative)
c 08/02/2016 V4.14 introduce upper limit for trace displays world
c coordinate range in counts; this makes the
c program robust, when displaying a zero offset
c trace with offset dependend scaling
c
c==============================================================================
c
......@@ -59,7 +66,7 @@ c
c
character*79 version
parameter(version=
& 'REFRACT V4.12 REFRACTion seismics - data interpretation')
& 'REFRACT V4.14 REFRACTion seismics - data interpretation')
c
c get common blocks
include 'refract_dim.inc'
......@@ -113,6 +120,7 @@ c read input data files
print *,'DEBUG: nfiles ',nfiles
endif
call setscale
call refract_setrefrange
c
call mpcfactors
call setfullrange
......
c this is <sub/refract_avgamp.f>
c ----------------------------------------------------------------------------
c
c Copyright (c) 2015 by Thomas Forbriger (BFO Schiltach)
c
c calculate average amplitude in offset range (used by subroutine mpcfactors)
c
c ----
c refract is free software; you can redistribute it and/or modify
c it under the terms of the GNU General Public License as published by
c the Free Software Foundation; either version 2 of the License, or
c (at your option) any later version.
c
c refract is distributed in the hope that it will be useful,
c but WITHOUT ANY WARRANTY; without even the implied warranty of
c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
c GNU General Public License for more details.
c
c You should have received a copy of the GNU General Public License
c along with this program. If not, see <http://www.gnu.org/licenses/>.
c ----
c
c REVISIONS and CHANGES
c 07/12/2015 V1.0 Thomas Forbriger
c
c ============================================================================
c
cS
real function refract_avgamp(ifile, allfiles)
c
integer ifile
logical allfiles
c
c Return average amplitude for traces in file ifile in amplitude
c reference range.
c
c ifile: file index of file to calculate average amplitude for
c allfiles: use all files not just file ifile
c
c return: average scaled amplitude value
cE
c
include 'refract_dim.inc'
include 'refract_data.inc'
include 'refract_para.inc'
c
integer navg, j
real maxamp, thisamp
logical usefile
c
navg=0
maxamp=0.
do j=1,ntraces
usefile=allfiles
if (.not.usefile) usefile=(fileindex(j).eq.ifile)
if (usefile
& .and.(fieldoffset(j).ge.plpar_avgrefxmin)
& .and.(fieldoffset(j).le.plpar_avgrefxmax)) then
navg=navg+1
if (plpar_remav) then
thisamp=max(abs(maxval(j)-average(j)),
& abs(minval(j)-average(j)))
else
thisamp=max(abs(maxval(j)),abs(minval(j)))
endif
maxamp=maxamp+thisamp*(fieldoffset(j)**plpar_expo)
endif
enddo
if (navg.lt.1) then
if (allfiles) then
print *,'ERROR (refract_avgamp): ',
& 'no traces in offset range selected by -S3 or -SN'
else
print *,'ERROR (refract_avgamp): ',
& 'file ',filename(ifile)(1:index(filename(ifile),' ')-1),
& ' has no traces in offset range selected by -S3 or -SN'
endif
stop 'aborting...'
endif
refract_avgamp=maxamp/navg
c
return
end
c
c ----- END OF sub/refract_avgamp.f -----
......@@ -39,6 +39,7 @@ c 24/10/2013 V1.12 - new option -So
c - new option -S3
c 18/11/2013 V1.13 new option -TF
c 21/03/2104 thof: new option -TR
c 07/12/2015 thof: new option -SN
c
c==============================================================================
c
......@@ -64,7 +65,7 @@ c declare local variables
integer i
c commandline
integer maxopt
parameter(maxopt=68)
parameter(maxopt=69)
character*3 optid(maxopt)
character*120 optarg(maxopt)
logical optset(maxopt), opthasarg(maxopt)
......@@ -136,7 +137,10 @@ c seismograms scaling
data optarg(65) /'0'/
data optid(66) /'-S3'/
data opthasarg(66) /.TRUE./
data optarg(66) /'0.,100.'/
data optarg(66) /'0.,10.'/
data optid(69) /'-SN'/
data opthasarg(69) /.TRUE./
data optarg(69) /'0.,10.'/
c
c additionals
data (optid(i), i=55,56) /'-Lt','-Ta'/
......@@ -221,8 +225,10 @@ c seismogram scaling
read(optarg(65), *) opt_Sordinate
if ((opt_Sordinate.lt.0).or.(opt_Sordinate.gt.3))
& stop 'ERROR: argument to -So is out of range'
opt_Savgref=optset(66)
read(optarg(66), *), opt_Savgrefxmin, opt_Savgrefxmax
opt_Savgref=optset(66).or.optset(69)
read(optarg(69), *), opt_Savgrefxmin, opt_Savgrefxmax
if (optset(66))
& read(optarg(66), *), opt_Savgrefxmin, opt_Savgrefxmax
if ((opt_Savgrefxmin.ge.opt_Savgrefxmax)
& .or.(opt_Savgrefxmin.lt.0.)
& .or.(opt_Savgrefxmax.lt.0.))
......@@ -261,11 +267,11 @@ c ============================
c
if (optset(24)) plflag_linestyle=opt_Lcycle
if (optset(29)) plflag_color=opt_Ccycle
if (optset(46)) plpar_expo=opt_Sexp
if (optset(46)) plpar_expo=max(0.,opt_Sexp)
plpar_forceminoff=optset(53)
if (plpar_forceminoff) plpar_minoff=opt_Sminoff
if (optset(54)) plpar_radius=opt_Sradius
if (opt_Savgref) opt_Smode=3
if (optset(66)) opt_Smode=3
plpar_mode=opt_Smode
print *,plpar_mode
c
......@@ -286,9 +292,8 @@ c
plflag_osnoreduce=opt_Sosnoreduce
plpar_vred=opt_Svel
c
plflag_m3avg=opt_Savgref
plpar_m3avgxmin=opt_Savgrefxmin
plpar_m3avgxmax=opt_Savgrefxmax
plpar_avgrefxmin=opt_Savgrefxmin
plpar_avgrefxmax=opt_Savgrefxmax
c
plflag_grid=opt_Egrid
plflag_vara=opt_Ewiggle
......
......@@ -381,7 +381,7 @@ c
flag_replot=.true.
elseif (plpar_pmmode.eq.3) then
if (verbose) print *,'reduce scaling exponent...'
plpar_expo=plpar_expo-0.08
plpar_expo=max(0.,plpar_expo-0.08)
call mpcfactors
flag_replot=.true.
else
......
......@@ -31,6 +31,8 @@ c offset for scaling mode 3 as the actual offset
c dependency is unknown
c 20/11/12 V1.2 use field offset for scaling purposes, not plot
c offset
c 07/12/15 V1.3 mode 2 and 3 now refer scaling to average scaled
c amplitude in defined offset range
c
c==============================================================================
cS
......@@ -45,22 +47,22 @@ c
include 'refract_para.inc'
c
cE
integer i,j,tref
real trefoff, refmpc, maxamp, thisamp
c real avgoff
integer i,j
real refmpc, refamp
real refract_avgamp
c
c scaling mode 1: individual scaling
c ----------------------------------
if (plpar_mode.eq.1) then
do j=1,ntraces
if (plpar_remav) then
maxamp=max(abs(maxval(j)-average(j)),abs(minval(j)-average(j)))
refamp=max(abs(maxval(j)-average(j)),abs(minval(j)-average(j)))
else
maxamp=max(abs(maxval(j)),abs(minval(j)))
refamp=max(abs(maxval(j)),abs(minval(j)))
endif
trv_mpc(j)=plpar_amp/maxamp
trv_mpc(j)=plpar_amp/refamp
if (debug) print *,'DEBUG (mpcfactors): trace, mpc ',j,trv_mpc(j)
if (debug) print *,'DEBUG (mpcfactors): maxamp ',maxamp
if (debug) print *,'DEBUG (mpcfactors): refamp ',refamp
if (debug) print *,'DEBUG (mpcfactors): maxval ',maxval(j)
if (debug) print *,'DEBUG (mpcfactors): minval ',minval(j)
enddo
......@@ -68,92 +70,22 @@ c
c scaling mode 2: least offset trace is reference
c -----------------------------------------------
elseif (plpar_mode.eq.2) then
c find trace with least fieldoffset
tref=1
trefoff=fieldoffset(1)
do j=1,ntraces
if (fieldoffset(j).lt.trefoff) then
tref=j
trefoff=fieldoffset(j)
endif
enddo
c calculate reference scale
if (plpar_remav) then
maxamp=max(abs(maxval(tref)-average(tref)),
& abs(minval(tref)-average(tref)))
else
maxamp=max(abs(maxval(tref)),abs(minval(tref)))
endif
refmpc=plpar_amp/(fieldoffset(tref)**plpar_expo)/maxamp
j=1
refamp=refract_avgamp(j, .true.)
refmpc=plpar_amp/refamp
c set mpc factors
do j=1,ntraces
trv_mpc(j)=refmpc*(fieldoffset(j)**plpar_expo)
enddo
c
c scaling mode 3: least offset trace in dataset is reference
c ----------------------------------------------------------
c scaling mode 3: least offset trace or average in given offset range
c in dataset is reference
c -------------------------------------------------------------------
elseif (plpar_mode.eq.3) then
c go for files
do i=1,nfiles
c print *,'ifile ',i
if (plflag_m3avg) then
c refer scaling to amplitude average
navg=0
maxamp=0.
c avgoff=0.
do j=1,ntraces
c print *,'itrace ',j,' offset ',fieldoffset(j)
if ((fileindex(j).eq.i)
& .and.(fieldoffset(j).ge.plpar_m3avgxmin)
& .and.(fieldoffset(j).le.plpar_m3avgxmax)) then
navg=navg+1
if (plpar_remav) then
thisamp=max(abs(maxval(j)-average(j)),
& abs(minval(j)-average(j)))
else
thisamp=max(abs(maxval(j)),abs(minval(j)))
endif
maxamp=maxamp+thisamp*(fieldoffset(j)**plpar_expo)
c avgoff=avgoff+fieldoffset(j)
endif
enddo
if (navg.lt.1) then
print *,'ERROR (mpc): ',
& 'file ',filename(i)(1:index(filename(i),' ')-1),
& ' has no traces in offset range selected by -S3'
stop 'aborting...'
endif
maxamp=maxamp/navg
c avgoff=avgoff/navg
refmpc=plpar_amp/maxamp
else
c refer scaling to nearest offset trace
c find trace with least fieldoffset within file
tref=0
do j=1,ntraces
if (fileindex(j).eq.i) then
if (tref.eq.0) then
tref=j
trefoff=fieldoffset(j)
else
if (fieldoffset(j).lt.trefoff) then
tref=j
trefoff=fieldoffset(j)
endif
endif
endif
enddo
if (debug) print *,'DEBUG: file ',i,' tref ',tref,
& ' trefoff ',trefoff
c calculate reference scale with respect to total reference
if (plpar_remav) then
maxamp=max(abs(maxval(tref)-average(tref)),
& abs(minval(tref)-average(tref)))
else
maxamp=max(abs(maxval(tref)),abs(minval(tref)))
endif
refmpc=plpar_amp/(maxamp*(fieldoffset(tref)**plpar_expo))
endif
refamp=refract_avgamp(i, .false.)
refmpc=plpar_amp/refamp
c set mpc factors within file
do j=1,ntraces
if (fileindex(j).eq.i) then
......
c this is <sub/refract_setrefrange.f>
c ----------------------------------------------------------------------------
c
c Copyright (c) 2015 by Thomas Forbriger (BFO Schiltach)
c
c Set offset range to take amplitude scaling reference from, if not set
c on command line.
c
c ----
c refract is free software; you can redistribute it and/or modify
c it under the terms of the GNU General Public License as published by
c the Free Software Foundation; either version 2 of the License, or
c (at your option) any later version.
c
c refract is distributed in the hope that it will be useful,
c but WITHOUT ANY WARRANTY; without even the implied warranty of
c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
c GNU General Public License for more details.
c
c You should have received a copy of the GNU General Public License
c along with this program. If not, see <http://www.gnu.org/licenses/>.
c ----
c
c REVISIONS and CHANGES
c 07/12/2015 V1.0 Thomas Forbriger
c
c ============================================================================
c
cS
subroutine refract_setrefrange
c
c set default scaling
c
include 'refract_dim.inc'
include 'refract_data.inc'
include 'refract_para.inc'
include 'refract_opt.inc'
cE
c
c The default is read from the command line parameter defaults. Make
c sure that the range at least includes on trace of each file.
c
real rmin,rmax
integer i
real reffraction
c fraction of total offset range to be used for reference
parameter(reffraction=0.1)
c
c adjust offset range only if not specified on command line
if (.not.opt_Savgref) then
c
c find overall minimum and maximum field offset value
rmin=fieldoffset(1)
rmax=fieldoffset(1)
do i=1,ntraces
rmin=min(rmin,fieldoffset(i))
rmax=max(rmax,fieldoffset(i))
enddo
c
c set reference to the smalles reffraction fraction
plpar_avgrefxmin=rmin
plpar_avgrefxmax=rmin+(rmax-rmin)*reffraction
c
c check for at least one trace of each file being present in the
c selected range
do i=1,nfiles
rmin=plpar_avgrefxmin
do j=1,ntraces
if (fileindex(j).eq.i) then
rmin=min(rmin,fieldoffset(j))
endif
enddo
plpar_avgrefxmax=max(plpar_avgrefxmax,rmin)
enddo
c
if (verbose) then
print *,'offset range to which amplitude scaling is referred:'
print *,' ',plpar_avgrefxmin,' m to ',plpar_avgrefxmax,' m'
endif
endif
c
return
end
c
c ----- END OF sub/refract_setrefrange.f -----
......@@ -5,6 +5,9 @@ c 30/04/98 by Thomas Forbriger (IfG Stuttgart)
c
c set default scaling
c
c This subroutine sets the initial amplitude, clip and minimum offset
c interval values.
c
c ----
c refract is free software; you can redistribute it and/or modify
c it under the terms of the GNU General Public License as published by
......
......@@ -27,6 +27,8 @@ c 03/07/98 V1.1 new clear concept - leave clipping to pgplot
c 04/07/98 V1.2 introduced traveltime reduction
c 12/11/12 V1.3 do not apply travel time reduction, if offset of
c seismograms is shifted intentionally
c 08/02/16 V1.4 introduce upper limit for world coordinate range
c in counts
c
c==============================================================================
cS
......@@ -37,10 +39,21 @@ c
c set viewport and world coordinate values for trace i
c
c returns true in case there is at least one sample to be plot
c
c reverse sign of scale if flag set
c
c values set are those required by calls
c call pgsvp(trv_vpleft, trv_vpright, trv_vpbot, trv_vptop)
c call pgswin(trv_tmin, trv_tmax, trv_vmin, trv_vmax)
c in subroutine pgtrace
c
c return values are passed through common block refract_trv which is
c defined in refract_seipar.inc
c
integer i
logical reverse
real limit, mpc, rmaxamp
parameter(limit=1.e38)
c
include 'refract_dim.inc'
include 'refract_data.inc'
......@@ -48,6 +61,22 @@ c
include 'refract_para.inc'
c
cE
c
c ----------------------------------------------------------------------
c limit
c -----
c
c To prevent numerical overflow, we must set an upper limit for floating
c point single precision variables. Fortran 95 provides a function
c huge(), where
c
c huge(1.0)=3.40282347E+38
c
c in a standard Intel environment with gfortran. Since this function
c might not be available for users relying on the code being compiled
c with a standard Fortran 77 compiler, I set the limit to 1.e38
c
c ----------------------------------------------------------------------
c
c tov stands for total viewport
c trv stands for trace viewport
......@@ -61,44 +90,151 @@ c
c the result flag is used to check whether any part of the curve is
c visible within the global viewport
c
c the baseline of the curve must be 0. or average
c the baseline of the curve must be 0. or averagec
c
c ----------------------------------------------------------------------
c
c The function calculates the viewport coordinates for the rectangle on
c the graphics device display to which the given trace signal shall be
c plotted, the trace viewport. These viewport coordinates are:
c
c left: trv_vpleft in normalized device coordinates
c right: trv_vpright in normalized device coordinates
c top: trv_vptop in normalized device coordinates
c bottom: trv_vpbot in normalized device coordinates
c
c The total viewport for the complete graphical display of traces is
c defined by the rectangle
c
c left: tov_vpleft in normalized device coordinates
c right: tov_vpright in normalized device coordinates
c top: tov_vptop in normalized device coordinates
c bottom: tov_vpbot in normalized device coordinates
c
c The subroutine make sure that the trace viewport is inside the total
c viewport.
c
c ----------------------------------------------------------------------
c
c The world coordinates for the total viewport are:
c
c largest offset: tov_rmax in meters
c smallest offset: tov_rmin in meters
c largest time: tov_tmax in seconds
c smallest time: tov_tmin in seconds
c
c With respect to these world coordinates in total viewport the world
c coordinates of the trace viewport are
c
c largest offset: trv_rtop in meters
c smallest offset: trv_rbot in meters
c largest time: trv_tright in seconds
c smallest time: trv_tleft in seconds
c
c ----------------------------------------------------------------------
c
c The actual world coordinates to be set to the window covering the
c trace viewport are:
c
c left: trv_tmin in counts
c right: trv_tmax in counts
c top: trv_vmax in seconds
c bottom: trv_vmin in seconds
c
c ----------------------------------------------------------------------
c
c The subroutine first sets the vertical extend and window in the offset
c domain and world coordinates in counts. In the second part it operates
c on horizontal coordinates and world coordinates in seconds.
c
c ----------------------------------------------------------------------
c
c
logical result
real vppm, ttmin, ttmax, theoffset
real realtime, plotoffset
c
c ======================================================================
c
c 1. Operate on vertical coordinates, offset domain and world
c coordinates in counts.
c -----------------------------------------------------------
c