Commit aaf8726e authored by thomas.forbriger's avatar thomas.forbriger Committed by thomas.forbriger
Browse files

reintegrated branch su1 on trunk

This is a legacy commit from before 2015-03-01.
It may be incomplete as well as inconsistent.
See COPYING.legacy and README.history for details.


SVN Path:     http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/trunk
SVN Revision: 3700
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parents 1fa56343 ee45f99f
......@@ -43,6 +43,9 @@
# hivgabor: calculate H/V from gabor output
# hivexpanco: claculate H/V from Fourier Bessel coefficients
#
# corresponding binaries being able to read SU and other formats:
# gaborx, phadix, gredax
#
# libraries/packages needed to compile the code
# liblapack: linear algebra package
# libblas: basic linear algebra functions
......@@ -54,11 +57,19 @@
# libts: package containing libts.a
# libsffu: package containing libsffu.a
# you can obtain the latter packages from where you obtained the present code
#
# to link binaries able to read SU data additional libraries are required:
# libfapidxx.a
# libdatrwxx.a
# libsffxx.a
# libgsexx.a
# libtime++.a
# libaff.a
#
# ============================================================================
#
all: greda gabor phadi hivexpanco hivgabor grereso
all: greda gabor phadi hivexpanco hivgabor grereso gaborx phadix gredax
#----------------------------------------------------------------------
# standard edit targets
......@@ -116,4 +127,19 @@ hivexpanco hivgabor gabor phadi: %: %.o
-L$(LOCLIBDIR) -L$(SERVERLIBDIR)
/bin/mv -fv $@ $(LOCBINDIR)
# binaries linked against libfapidxx and colleagues, being able to read SU
# data and other formats:
gredax: greda.o greda_phasor.o
$(FC) -o $@ $^ -lsffu -lts -ltf $(LINLIB) \
-lgsl -lgslcblas -ltime -lfapidxx -ldatrwxx -lsffxx -lgsexx\
-ltime++ -laff \
-L$(LOCLIBDIR) -L$(SERVERLIBDIR)
/bin/mv -fv $@ $(LOCBINDIR)
gaborx phadix: %x: %.o
$(FC) -o $@ $< -lsffu -lts -ltf $(LINLIB) \
-ltime -lfapidxx -ldatrwxx -lsffxx -lgsexx -laff -ltime++ \
$(PGPLOTLIB) \
-L$(LOCLIBDIR) -L$(SERVERLIBDIR)
/bin/mv -fv $@ $(LOCBINDIR)
# ----- END OF Makefile -----
......@@ -27,22 +27,23 @@ c 22/02/2001 V1.0 Thomas Forbriger
c 02/08/2001 V1.1 allow trace order rather than offset order
c 09/09/2004 V1.2 set lower limit of time window
c 13/06/2006 V1.3 provide Hanning taper
c 30.12.2010 V1.4 implemented file format selection
c
c==============================================================================
c
program gabor
c
character*(*) version
parameter(version='GABOR V1.3 calculate gabor matrix')
parameter(version='GABOR V1.4 calculate gabor matrix')
character*(*) GABOR_CVS_ID
parameter(GABOR_CVS_ID=
& '$Id$')
c
c input dataset
character*80 filename
character*80 filename, informat
integer maxtraces, totmaxsamples
parameter(maxtraces=40, totmaxsamples=4000000)
integer lu
integer lu, ierr
parameter(lu=12)
real fdata(totmaxsamples)
integer idata(totmaxsamples)
......@@ -87,15 +88,17 @@ c debugging
c commandline
integer maxopt, lastarg, iargc
character*80 argument
parameter(maxopt=11)
character*2 optid(maxopt)
parameter(maxopt=12)
character*3 optid(maxopt)
character*40 optarg(maxopt)
logical optset(maxopt), opthasarg(maxopt)
c here are the keys to our commandline options
data optid/2h-d, 2h-v, '-f', '-n', '-o', '-t', '-T', '-w', '-S',
& 2h-m, '-h'/
data opthasarg/2*.FALSE.,2*.TRUE.,.FALSE.,5*.TRUE.,.false./
data optarg/2*1h-,'100.','100','-','1.','1','1.','1','0.','-'/
data optid/'-d', '-v', '-f', '-n', '-o', '-t', '-T', '-w', '-S',
& '-m', '-h', '-ty'/
data opthasarg/2*.FALSE.,2*.TRUE.,.FALSE.,5*.TRUE.,.false.,
& .true./
data optarg/2*'-','100.','100','-','1.','1','1.','1','0.','-',
& 'sff'/
c
c------------------------------------------------------------------------------
c basic information
......@@ -103,14 +106,20 @@ c
c
argument=' '
if (iargc().eq.1) call getarg(1, argument)
c
if (argument(1:6).eq.'-xhelp') then
call sff_help_details
stop
endif
c print *,iargc()
if ((argument(1:5).eq.'-help').or.(iargc().lt.2)) then
print *,version
print *,'Usage: gabor datafile outfile [-v] [-d]'
print *,' [-f freq] [-n N] [-t tmax] [-o]'
print *,' [-T N | -S N] [-h]'
print *,' [-w fac] [-m tmin]'
print *,' [-w fac] [-m tmin] [-ty f]'
print *,' or: gabor -help'
print *,' or: gabor -xhelp'
if (argument(1:5).ne.'-help')
& stop 'ERROR: wrong number of arguments'
print *,' '
......@@ -132,8 +141,11 @@ c print *,iargc()
print *,'-w fac broaden time windoe by factor ''fac'' '
print *,'-m tmin use windows starting at tmin'
print *,'-h use Hanning taper rather than Gaussian'
print *,'-ty f select input file format'
print *,' '
print *,GABOR_CVS_ID
print *,' '
call sff_help_formats
stop
endif
c
......@@ -154,6 +166,7 @@ c
if (optset(9)) read(optarg(9), *) nptraces
read(optarg(10), *) tmin
hanning=optset(11)
informat=optarg(12)
overwrite=optset(5)
......@@ -163,6 +176,8 @@ c
call getarg(2, outbase)
c------------------------------------------------------------------------------
c go
call sff_select_input_format(informat, ierr)
if (ierr.ne.0) stop 'ERROR: selecting input file format'
call sffu_simpleread(lu, filename, maxtraces, totmaxsamples,
& fdata, idata, toffset, tracedt, roffset, innsamples,
& firstsample, ntraces, verbose)
......
......@@ -103,7 +103,8 @@ c spectra is sample index
c 13/11/10 V3.20 - do not expose misleading term "green" to the
c user
c - use GSL instead of numerical recipes
c 12/01/2011 V3.20b - added definition of Fourier transformation
c 30.12.2010 V3.21 - implemented libfapidxx
c 12/01/2011 V3.21b - added definition of Fourier transformation
c online help text
c
c==============================================================================
......@@ -114,7 +115,7 @@ c first we declare some general variables
c
character*79 version,CVSID
parameter(version=
&'GREDA V3.20b Calculate Fourier-Bessel expansion coefficients')
&'GREDA V3.21b Calculate Fourier-Bessel expansion coefficients')
parameter(CVSID=
& '$Id$')
c
......@@ -130,7 +131,7 @@ c----------------------------------------------------------------------
c here follows everything we need for the input data hold
c
c datafile
character*80 filename, Fourierfile
character*80 filename, Fourierfile, informat
real fdata(maxsamp)
integer idata(maxsamp)
equivalence(fdata, idata)
......@@ -186,27 +187,28 @@ c parameters
logical pwofile, pwoautofile, pwosel
c command line
integer maxopt, lastarg, iargc
parameter(maxopt=35)
parameter(maxopt=36)
character*4 optid(maxopt)
character*80 optarg(maxopt)
logical optset(maxopt), opthasarg(maxopt)
c here are the keys to our commandline options
data optid/2h-d,2h-t,2h-s,2h-n,2h-f,2h-o,2h-1,2h-2,2h-v,2h-T,2h-S,
& 2h-M,2h-L,2h-q,2h-E,2h-K,2h-D,2h-a,2h-b,2h-g,2h-H,2h-O,2h-Q,2h-B,
& 2h-P,2h-r,2h-F,2h-N,2h-X,2h-W,4h-tap,4h-spo,3h-pw,4h-pwf,4h-pwa/
& 2h-P,2h-r,2h-F,2h-N,2h-X,2h-W,4h-tap,4h-spo,3h-pw,4h-pwf,4h-pwa,
& '-ty'/
data opthasarg/.FALSE.,4*.TRUE.,4*.FALSE.,.TRUE.,3*.FALSE.,2*.TRUE.,
& 6*.FALSE.,2*.TRUE.,.TRUE.,.FALSE.,.TRUE.,.FALSE.,.TRUE.,.FALSE.,
& 6*.TRUE./
& 7*.TRUE./
data optarg/1h-,3h10.,2h8.,1h5,4h100.,4*1h-,2h0.,3*1h-,5h1.,1.,
& 2h1.,6*1h-,4h0.01,5h1.,1.,4*1h-,3h10.,1h-,6h1.,10.,
& 11h0.,0.,0.,0.,4*3hnil/
& 11h0.,0.,0.,0.,4*3hnil,'sff'/
c
c======================================================================
c
c go on with executable code
c
print *,version
print *,'Usage: greda datafile coeffile'
print *,'Usage: greda datafile coeffile [-ty format]'
print *,' [-L] [-K] [-P] [-D] [-M] [-H]'
print *,' [-n nslo] [-s smax] [-f fmax]'
print *,' [-t frac] [-T frac] [-E edge]'
......@@ -220,10 +222,14 @@ c
print *,' [-pwf filename]'
print *,' [-pwa filename]'
print *,'or greda -help'
print *,'or greda -xhelp'
c
if (iargc().lt.1) stop 'ERROR: missing parameters'
call getarg(1, filename)
if (filename(1:5).eq.'-help') then
if (filename(1:6).eq.'-xhelp') then
call sff_help_details
stop
else if (filename(1:5).eq.'-help') then
print *,' '
print *,'Calculate Fourier-Bessel expansion coefficients'
print *,'for a single-shot seismic profile'
......@@ -239,6 +245,7 @@ c
print *,' the same time of first sample and the same'
print *,' sampling interval.'
print *,'coeffile Name of file to contain results.'
print *,'-ty format input file format (see list below)'
print *,' '
print *,'Available alternatives to the Fourier-Bessel transform:'
print *,'-L Calculate by linear inversion (exp-damping).'
......@@ -425,6 +432,8 @@ c
print *,'magic numbers:'
print *,' Fourier-Bessel coefficient files: ',cmagic
print *,' Trace coefficient files: ',cspmagic
print *,' '
call sff_help_formats
stop
endif
c
......@@ -514,6 +523,7 @@ c print *,'DEBUG: returned from tf_cmdline'
pwffilename=optarg(34)
pwoautofile=optset(35)
pwafilename=optarg(35)
informat=optarg(36)
if (debug) print *,'DEBUG: read options'
c
if ((.not.(planewave)).and.(smax.lt.0.))
......@@ -558,6 +568,8 @@ c go for the real calculations
c
c read seismic data
if (debug) print *,'DEBUG: call readdata'
call sff_select_input_format(informat, ierr)
if (ierr.ne.0) stop 'ERROR: selecting input file format'
call readdata(filename, fdata, idata, spectra, r, maxr,
& maxtr, maxsamp, ntr, nsamp, dt, tfirst)
if (debug) print *,'DEBUG: returned from readdata'
......@@ -1112,7 +1124,7 @@ c open sff-file
if ((ntr.lt.maxtr).and.(.not.last)) goto 1
if (.not.(last)) then
print *,'NOTICE: will ignore traces after ',ntr
close(lu)
call sff_close(lu)
endif
print *,'file read and closed'
return
......
......@@ -91,14 +91,16 @@ c
if (argument(1:5).ne.'-help')
& stop 'ERROR: wrong number of arguments'
print *,' '
print *,'Hinfile horizontal component input file'
print *,'Vinfile vertical component input file'
print *,'Hinfile horizontal component input file from syg'
print *,'Vinfile vertical component input file from syg'
print *,'outfile output file name'
print *,' '
print *,'The program calculates the complex H/V ratio of for'
print *,'Fourier Bessel coefficients as computed by sug.'
print *,'All files are in the format produced by syg.f'
print *,' '
print *,'Results can be displayed by grepg'
print *,' '
print *,'-v be verbose'
print *,'-o overwrite existing output file'
print *,'-w val water level for denominator as a fraction'
......
......@@ -90,9 +90,9 @@ c
if (argument(1:5).ne.'-help')
& stop 'ERROR: wrong number of arguments'
print *,' '
print *,'Hinfile horizontal component input file'
print *,'Vinfile vertical component input file'
print *,'outfile output file name'
print *,'Hinfile horizontal component spectrogram file'
print *,'Vinfile vertical component spectrogram file'
print *,'outfile output spectrogram file name'
print *,' '
print *,'The program calculates the complex ratio of gabor'
print *,'H/V from a horizontal component and vertical component'
......
......@@ -34,6 +34,7 @@ c 25/07/2010 V1.5 Zacharias found a bug: There was a parameter
c missing in a call to pgpt - I wonder how this
c went unnoticed over so many years...
c 27/07/2010 V1.6 correct plot labels
c 30.12.2010 V1.7 implemented libfapidxx
c
c==============================================================================
c
......@@ -41,15 +42,15 @@ c
c
character*(*) version
parameter(version=
& 'PHADI V1.6 calculate dispersion relation from phase differences')
& 'PHADI V1.7 calculate dispersion relation from phase differences')
character*(*) PHADI_CVS_ID
parameter(PHADI_CVS_ID=
& '$Id$')
c input dataset
character*80 filename
character*80 filename, informat
integer maxtraces, totmaxsamples
parameter(maxtraces=300, totmaxsamples=2000000)
integer lu
integer lu, ierr
parameter(lu=12)
real data(totmaxsamples)
integer idata(totmaxsamples)
......@@ -98,7 +99,7 @@ c constants
c commandline
integer maxopt, lastarg, iargc
character*80 argument
parameter(maxopt=22)
parameter(maxopt=23)
character*3 optid(maxopt)
character*80 optarg(maxopt)
logical optset(maxopt), opthasarg(maxopt)
......@@ -107,12 +108,12 @@ c debugging
c here are the keys to our commandline options
data optid/3h-DE, 2h-v, '-P', '-Pf', '-Po', '-d', '-p', '-pf', '-po',
& '-pt', '-D', '-Df', '-Do', '-Dv', '-f', '-Dt',
& '-DO', '-Dn', '-Dr', '-Dp','-De', '-Dl'/
& '-DO', '-Dn', '-Dr', '-Dp','-De', '-Dl', '-ty'/
data opthasarg/3*.FALSE.,3*.TRUE.,.FALSE.,2*.TRUE.,2*.FALSE.,
& 2*.TRUE.,.FALSE.,5*.TRUE.,2*.FALSE.,.true./
& 2*.TRUE.,.FALSE.,5*.TRUE.,2*.FALSE.,2*.true./
data optarg/3*1h-,'10.,20.','phasor.out','x11','-','10.,20.',
& 'phase.out',2*'-','10.,20.','disan.out','-','10.,20.',
& '0.', '1.e-4', '10', '1,1',2*'-','0.'/
& '0.', '1.e-4', '10', '1,1',2*'-','0.','sff'/
c
c------------------------------------------------------------------------------
c basic information
......@@ -128,7 +129,13 @@ c
print *,' [-D] [-Df f1,f2] [-Do file] [-Dv]'
print *,' [-Dt val] [-DO eps] [-Dn N]'
print *,' [-Dr Nl,Nr] [-Dp] [-De] [-Dl f]'
print *,' [-ty format]'
print *,' or: phadi -help'
print *,' or: phadi -xhelp'
if (argument(1:6).eq.'-xhelp') then
call sff_help_details
stop
endif
if (argument(1:5).ne.'-help') stop 'ERROR: wrong number of arguments'
print *,' '
print *,'calculate dispersion relation from phase differences'
......@@ -139,6 +146,7 @@ c
print *,'-DE print debuggin output'
print *,'-d device PGPLOT output device'
print *,'-f f1,f2 default frequency interval'
print *,'-ty format select file format (see list below)'
print *,' '
print *,'-P execute phasor test'
print *,'-Pf f1,f2 frequency interval for phasor test'
......@@ -172,6 +180,8 @@ c
print *,' '
call pgp_showdevices
print *,' '
call sff_help_formats
print *,' '
print *,PHADI_CVS_ID
stop
endif
......@@ -215,11 +225,14 @@ c
fitplot=optset(20)
errorplot=optset(21)
read(optarg(22), *) disanthreshlim
informat=optarg(23)
call getarg(1, filename)
c
c------------------------------------------------------------------------------
c go
call sff_select_input_format(informat, ierr)
if (ierr.ne.0) stop 'ERROR: selecting input file format'
call sffu_simpleread(lu, filename, maxtraces, totmaxsamples, data, idata,
& toffset, tracedt, roffset, innsamples, firstsample, ntraces, verbose)
c
......
......@@ -29,11 +29,12 @@
*
* REVISIONS and CHANGES
* - 19/11/2010 V1.0 Thomas Forbriger
* - 23.12.2010 V1.1 use stream buffer seekg() function to skip series
*
* ============================================================================
*/
#define DATRW_SU_CC_VERSION \
"DATRW_SU_CC V1.0 "
"DATRW_SU_CC V1.1"
#define DATRW_SU_CC_CVSID \
"$Id$"
......@@ -92,7 +93,12 @@ namespace datrw {
void isustream::skipseries()
{
Tfseries s=this->fseries();
Mis.seekg(sizeof(float)*Mnextheader.Mheader.ns,
std::ios_base::cur);
DATRW_Xassert(Mis.good(),
"ERROR (isustream::skipseries): calling seekg()",
::datrw::su::SUReadException);
this->readheader();
} // void isustream::skipseries()
/*----------------------------------------------------------------------*/
......
......@@ -30,6 +30,19 @@
*
* REVISIONS and CHANGES
* - 06/12/2010 V1.0 Thomas Forbriger
* - 14/01/2011 V1.1 extended modifications:
* - resolved problem with coordinate values close to
* zero which lead to a logarithm approaching
* infinity upon testing the number of significant
* digits
* - resolved a problem with unreasonable large scale
* values in case of small round off errors in
* floating point coordinate values
* - apply proper rounding to coordinate values, such
* that 0.0999999 will be stored as 0.1
* - check actual round-off error rather than scale
* and issue a warning in cases where round-off error
* exceed 0.1%
*
* ============================================================================
*/
......@@ -68,6 +81,11 @@ namespace datrw {
/*======================================================================*/
double ScalCoo::effectivezero=1.e-12;
int ScalCoo::maxnsigdigits=4;
/*----------------------------------------------------------------------*/
//! set from header fields
void ScalCoo::set(const short& s, const int& c)
{
......@@ -97,42 +115,80 @@ namespace datrw {
{
const bool debug=false;
this->scale=1;
this->coo=static_cast<int>(value);
double v=value>0 ? value : -value;
DATRW_debug(debug, "ScalCoo::set(const double& v)",
"v: " << v << "\n"
"log10(v): " << std::log10(v) << "\n"
"floor(log10(v)): " << std::floor(std::log10(v)));
int nlead=1+static_cast<int>(std::floor(std::log10(v)));
int ntrail=::datrw::util::ntrailingdigits(v);
int nsig=::datrw::util::nsignificantdigits(v, debug);
DATRW_debug(debug, "ScalCoo::set(const double& v)",
"nlead: " << nlead << "\n"
"ntrail: " << ntrail << "\n"
"nsig: " << nsig);
if (ntrail>0)
this->coo=static_cast<int>(nearbyint(value));
double v=static_cast<float>(value>=0 ? value : -value);
if (v<ScalCoo::effectivezero)
{
int s=std::pow(10,ntrail);
DATRW_report_assert(s<100000,
"ERROR (ScalCoo::set): scale too large",
"scale: " << s << " coordinate value "
<< helper::MyOutputFormat() << value
<< " will be truncated!")
s = s > 10000 ? 10000 : s;
DATRW_assert(s<100000,
"ERROR (ScalCoo::set): scale too large");
this->scale=-s;
this->coo=static_cast<int>(value*s);
this->scale=1;
this->coo=0;
}
else if (nlead>nsig)
else
{
int s=std::pow(10,nlead-nsig);
s = s > 10000 ? 10000 : s;
DATRW_assert(s<100000,
"ERROR (ScalCoo::set): scale too large");
this->scale=s;
this->coo=static_cast<int>(value/s);
}
DATRW_debug(debug, "ScalCoo::set(const double& v)",
"v: " << helper::MyOutputFormat() << v << "\n"
"log10(v): " << std::log10(v) << "\n"
"floor(log10(v)): " << std::floor(std::log10(v)));
int nlead=1+static_cast<int>(std::floor(std::log10(v)));
int ntrail=::datrw::util::ntrailingdigits(v);
int nsig=::datrw::util::nsignificantdigits(v, debug);
int nzero=ntrail-nsig;
DATRW_debug(debug, "ScalCoo::set(const double& v)",
"nlead: " << nlead << "\n"
"ntrail: " << ntrail << "\n"
"nsig: " << nsig);
nsig = nsig>maxnsigdigits ? maxnsigdigits : nsig;
ntrail=nsig+nzero;
DATRW_debug(debug, "ScalCoo::set(const double& v)",
"nlead: " << nlead << "\n"
"ntrail: " << ntrail << "\n"
"nsig: " << nsig);
if (ntrail>0)
{
DATRW_debug(debug, "ScalCoo::set(const double& value)",
"ntrail>0");
int s=std::pow(10,ntrail);
s = s > 10000 ? 10000 : s;
DATRW_assert(s<100000,
"ERROR (ScalCoo::set): scale too large");
this->scale=-s;
this->coo=static_cast<int>(nearbyint(value*s));
// check truncation error
double truevalue=value*s;
double roundedvalue=nearbyint(truevalue);
double truncationerror=fabs(1.-(truevalue/roundedvalue));
DATRW_debug(debug, "ScalCoo::set(const double& value)",
"truncation error:"
<< " truevalue: " << truevalue
<< " roundedvalue: " << roundedvalue
<< " truncationerror: " << truncationerror);
DATRW_report_assert(truncationerror<0.001,
"WARNING (ScalCoo::set) "
"truncation error > 0.1%",
"coordinate value "
<< helper::MyOutputFormat() << value
<< " will be truncated!"
<< "\n*"
<< " scaled value: " << truevalue
<< " rounded scaled value: "
<< roundedvalue);
}
else if (nlead>nsig)
{
DATRW_debug(debug, "ScalCoo::set(const double& value)",
"nlead>nsig");
int s=std::pow(10,nlead-nsig);
s = s > 10000 ? 10000 : s;
DATRW_assert(s<100000,
"ERROR (ScalCoo::set): scale too large");
this->scale=s;
this->coo=static_cast<int>(nearbyint(value/s));
}
} // if (v<ScalCoo::effectivezero), else clause
DATRW_debug(debug, "ScalCoo::set(const double& value)",
"value on exit:"
<< " value: " << value
<< " this->coo: " << this->coo
<< " this->scale: " << this->scale);
} // void ScalCoo::set(const double& v)
/*----------------------------------------------------------------------*/
......@@ -175,6 +231,8 @@ namespace datrw {
int pp = p < 0 ? -p : p;
short s=static_cast<short>(std::pow(10,pp));
this->scale= p<0 ? -s: s;
DATRW_assert(this-scale != 0,
"ERROR (ScalCoo::scaletopower): illegal scale value!");
} // void ScalCoo::scaletopower(const int& p)
/*----------------------------------------------------------------------*/
......@@ -196,6 +254,8 @@ namespace datrw {
int pp = p < 0 ? -p : p;
short s=static_cast<short>(std::pow(10,pp));
this->scale= p<0 ? -s: s;
DATRW_assert(this-scale != 0,
"ERROR (ScalCoo::scaletopower): illegal scale value!");
} // void ScalCoo::adjustscaling()
/*----------------------------------------------------------------------*/
......@@ -213,6 +273,17 @@ namespace datrw {
//! read values from SU header
void Coordinates::getvaluesfrom(const TraceHeaderStruct& h)
{