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

preovide Fourier input and output for phasor walkout

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.
greda: new spatial taper function


SVN Path:     http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/trunk
SVN Revision: 924
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 19859ed7
c this is <greda.f>
c------------------------------------------------------------------------------
c $Id: greda.f,v 1.11 2002-06-13 11:16:50 forbrig Exp $
c $Id: greda.f,v 1.12 2002-09-13 21:05:01 forbrig Exp $
c
c 24/06/97 by Thomas Forbriger (IfG Stuttgart)
c
......@@ -74,6 +74,8 @@ c planestack for the Slant Stack algorithm and into
c subroutine expmodel. In expmodel it makes only sense
c together with the Bessel transform, since only in that
c case the inverse of the Gram matrix is diagonal.
c 13/09/02 V3.13 - add extra offset taper
c - output seismogram spectra
c
c==============================================================================
c
......@@ -82,7 +84,7 @@ c
c first we declare some general variables
c
character*79 version
parameter(version='GREDA V3.12 Greens function from data')
parameter(version='GREDA V3.13 Greens function from data')
c
c calculations common block
include 'greda.inc'
......@@ -95,7 +97,7 @@ c----------------------------------------------------------------------
c here follows everything we need for the input data hold
c
c datafile
character*80 filename
character*80 filename, fourierfile
real fdata(maxsamp)
integer idata(maxsamp)
equivalence(fdata, idata)
......@@ -116,6 +118,10 @@ c magic number for binary file identification
integer magic
character*4 cmagic
parameter(cmagic='1234')
c magic number for trace spectra binary file identification
integer spmagic
character*4 cspmagic
parameter(cspmagic='SP34')
c greens function
character*80 greensfile
complex green(maxslo, maxom)
......@@ -131,27 +137,29 @@ c
c parameters
real tapfrac, offtapfrac, edgefrac, minoff, stackdelta, rescaleexpo
real hopnumberfrequency, wltaplen, wltapfrac
real tapoffsets(4)
logical hopnumberthis, backtranscale
logical debug, overwrite, hankel1, hankel2, verbose, uzerospecial
logical matrixmethod, lininv, offtaper, edgeset, linkinv
logical disgram, softcosine, gausstaper, gausstime, parkermethod
logical stackso, planewave, rescale, specrescale, hopnumbers
logical applywltaper
logical applywltaper,specialtaper,writefourier
c command line
integer maxopt, lastarg, iargc
parameter(maxopt=30)
character*2 optid(maxopt)
character*40 optarg(maxopt)
parameter(maxopt=32)
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/
& 2h-P,2h-r,2h-F,2h-N,2h-X,2h-W,4h-tap,4h-spo/
data opthasarg/.FALSE.,4*.TRUE.,4*.FALSE.,.TRUE.,3*.FALSE.,2*.TRUE.,
& 6*.FALSE.,2*.TRUE.,.TRUE.,.FALSE.,.TRUE.,.FALSE.,.TRUE.,.FALSE.,
& .TRUE./
& 3*.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./
& 2h1.,6*1h-,4h0.01,5h1.,1.,4*1h-,3h10.,1h-,6h1.,10.,
& 11h0.,0.,0.,0.,3hnil/
c
c======================================================================
c
......@@ -166,6 +174,8 @@ c
print *,' [-O minoff] [-B delta] [-r expo] [-F]'
print *,' [-1] [-2] [-q f,e] [-Q f,e] [-S]'
print *,' [-v] [-o] [-N f] [-X]'
print *,' [-tap o1,o2,o3,o4]'
print *,' [-spo filename]'
print *,'or greda -help'
c
if (iargc().lt.1) stop 'ERROR: missing parameters'
......@@ -174,7 +184,7 @@ c
print *,' '
print *,'Calculate Green''s Function from a set of seismograms'
print *,'by using linear inversion theory. Actually does'
print *,'a fourier-bessel back transform by default.'
print *,'a fourier-bessel transform by default.'
print *,' '
print *,'datafile Name of file containing seismograms.'
print *,' The program will expect a homogeneous dataset.'
......@@ -183,7 +193,7 @@ c
print *,' sampling interval.'
print *,'greensfile Name of file to contain results.'
print *,' '
print *,'Available alternatives to the fourier-bessel back transform:'
print *,'Available alternatives to the fourier-bessel transform:'
print *,'-L Calculate by linear inversion (exp-damping).'
print *,'-K Calculate by linear inversion (K_0-damping).'
print *,'-H Use the method following Henry, Orcutt and Parker:'
......@@ -227,6 +237,15 @@ c
print *,' compared to it''s original value.'
print *,' The default is not to apply any offset domain'
print *,' taper.'
print *,'-tap o1,o2,o3,o4'
print *,' special offset taper'
print *,' the taper is 0. for offsets smaller than'
print *,' o1 and offsets larger than o4'
print *,' the taper i 1. for offsets between o2'
print *,' and o3'
print *,' it has a sine edge between o1 and'
print *,' and o2 and a cosine edge between o3'
print *,' and o4'
print *,'-E edge Specify the sample from which on all'
print *,' samples should be tapered to zero as a'
print *,' fraction of the time series length.'
......@@ -293,6 +312,9 @@ c
print *,'-o Overwrite existing output file.'
print *,'-v Be somehow verbose.'
print *,'-N f print HOP numbers for frequency f'
print *,'-spo filename'
print *,' write seismic traces'' Fourier'
print *,' coefficients to file ''filename'' '
print *,' '
print *,'How it works:'
print *,'In cases -L, -K, -H and -D (linear inversion) we use'
......@@ -320,6 +342,10 @@ c
print *,' traces: ',maxtr
print *,' slownesses: ',maxslo
print *,' frequencies: ',maxom
print *,' '
print *,'magic numbers:'
print *,' Green coefficient files: ',cmagic
print *,' Trace coefficient files: ',cspmagic
stop
endif
c
......@@ -388,6 +414,11 @@ c
backtranscale=optset(29)
applywltaper=optset(30)
read(optarg(30), *) wltaplen, wltapfrac
specialtaper=optset(31)
if (specialtaper) read(optarg(31), *) (tapoffsets(i), i=1,4)
writefourier=optset(32)
fourierfile=optarg(32)
c
if ((.not.(planewave)).and.(smax.lt.0.))
& stop 'ERROR: negative maximum slowness'
c report
......@@ -444,6 +475,11 @@ c apply tapers
call taper(maxsamp, maxtr, ntr, nsamp, tapfrac, offtapfrac,
& maxr, spectra, r, offtaper, verbose, edgeset, edgefrac,
& softcosine, gausstaper, gausstime)
c
c apply special tapers
if (specialtaper)
& call specialtap(maxsamp, maxtr, ntr, nsamp, tapoffsets, maxr,
& spectra, r, verbose)
c
c rescale traces if we want them to be rescaled
if (rescale) call dorescale(rescaleexpo, r)
......@@ -451,6 +487,11 @@ c
c calculate complex spectra
call calcspec(dt, tfirst, om,
& fmax, optset(5))
c
c write Fourier coefficients
if (writefourier)
& call fourierwrite(fourierfile, spmagic, cspmagic, overwrite,
& r, maxr, om)
c
c rescale traces if we want them to be rescaled
if (specrescale) call dospecrescale(rescaleexpo, r, om)
......@@ -734,11 +775,24 @@ c subroutine taper(maxsamp, maxtr, ntr, nsamp, tapfrac, offtapfrac,
c & maxr, spectra, r, offtaper, verbose, edgeset, edgefrac,
c & softcosine, gausstaper, gausstime)
c
c apply special offset domain taper (option -tap tapoffsets)
c ---------------------------------
c
c subroutine specialtap(maxsamp, maxtr, ntr, nsamp, tapoffsets, maxr
c & spectra, r, verbose)
c
c calculate complex Fourier coefficients from waveform data
c ---------------------------------------------------------
c reads common blocks from greda.inc
c
c subroutine calcspec(dt, tfirst, om, fmax, fmaxset)
c
c write trace Fourier spectra (easy to use) (option -spo filename)
c ----------------------------------------
c reads common blocks from greda.inc
c
c subroutine fourierwrite(filename, magic, cmagic, overwrite,
c & r, maxr, om)
c
c set slowness values for different expansion cases
c -------------------------------------------------
......@@ -1172,6 +1226,56 @@ c
end
c
c----------------------------------------------------------------------
c
subroutine specialtap(maxsamp, maxtr, ntr, nsamp, tapoffsets,
& maxr, spectra, r, verbose)
c
c apply special offset domain taper
c
integer maxsamp, maxtr, ntr, nsamp
real tapoffsets(4), maxr
double complex spectra(maxtr, maxsamp)
real r(maxtr)
logical verbose
c
integer i, l
real fac, pi
parameter(pi=3.141592653589793)
c
print *,' '
print *,'apply special offset taper to input data'
print *,' min=',tapoffsets(1),', left=',tapoffsets(2),
& ' right=',tapoffsets(3),', max=',tapoffsets(4)
c
c apply offset domain taper
c
do i=1,ntr
fac=1.
if (r(i).lt.tapoffsets(1)) then
fac=0.
elseif (r(i).gt.tapoffsets(4)) then
fac=0.
elseif (r(i).lt.tapoffsets(2)) then
fac=sin(0.5*pi*(r(i)-tapoffsets(1))/
& (tapoffsets(2)-tapoffsets(1)))
elseif (r(i).gt.tapoffsets(3)) then
fac=cos(0.5*pi*(r(i)-tapoffsets(3))/
& (tapoffsets(4)-tapoffsets(3)))
else
fac=1.
endif
if (verbose) then
print *,' at ',r(i),'m fac: ',fac
endif
do l=1,nsamp
spectra(l,i)=fac*spectra(l,i)
enddo
enddo
c
return
end
c
c----------------------------------------------------------------------
c
subroutine calcspec(dt, tfirst, om, fmax, fmaxset)
c
......@@ -1259,6 +1363,58 @@ c
end
c
c----------------------------------------------------------------------
c
c write trace Fourier spectra (easy to use)
c
subroutine fourierwrite(filename, magic, cmagic, overwrite,
& r, maxr, om)
c
include 'greda.inc'
c
character*4 cmagic
integer magic
character*(*) filename
logical overwrite
real r(maxtr), maxr, om(maxom)
c
integer chain(maxtr), firstinchain,i,j
integer k(maxtr)
integer lu
parameter(lu=20)
c
print *,' '
if (overwrite) then
print *,'opening Fourier file ',filename(1:index(filename,' ')),
& ' - overwrite mode'
open(lu, file=filename, form='unformatted', err=98)
else
print *,'opening Fourier file ',filename(1:index(filename,' '))
open(lu, file=filename, status='new',
& form='unformatted', err=98)
endif
c
c sort offsets
call tf_rchain(r, chain, ntr, firstinchain, 1)
j=firstinchain
do i=1,ntr
k(i)=j
j=chain(j)
enddo
c
call tf_magic(cmagic, magic)
write(lu, err=97) magic
write(lu, err=97) nom, ntr
write(lu, err=97) (om(i), i=1,nom), (r(k(i)), i=1,ntr)
write(lu, err=97) ((spectra(k(j),i), i=1,nom), j=1,ntr)
close(lu, err=96)
c
return
98 stop 'ERROR: opening Fourier file'
97 stop 'ERROR: writing Fourier file'
96 stop 'ERROR: closing Fourier file'
end
c
c----------------------------------------------------------------------
c
subroutine setslo(maxslo, nslo, smax, slo,
& hankel1, hankel2, uzerospecial)
......
......@@ -13,7 +13,7 @@ CFLAGS=-O2
GREBOBS=grepg.o grepg_message.o grepg_dopicks.o grepg_selstyle.o \
grepg_phase.o grepg_phasewedg.o grepg_poly.o grepg_remavg.o \
grepg_contr.o grepg_prepcol.o
grepg_contr.o grepg_prepcol.o grepg_readfourier.o
make.incdep: *.f
incdep > make.incdep
......
c this is file <grepg.f>
c======================================================================
c $Id: grepg.f,v 1.25 2002-06-12 20:28:01 forbrig Exp $
c $Id: grepg.f,v 1.26 2002-09-13 21:05:02 forbrig Exp $
c
c GREPG.F
c
......@@ -87,6 +87,7 @@ c V2.33 22/01/02 allow color and linestyle cycle switching
c V2.33a 02/04/02 give more specific error messages
c V2.34 19/05/02 allow plots matching different Fourier transform sign
c conventions
c V2.35 13/09/02 phasor walkout from Fourier data
c
c======================================================================
program grepg
......@@ -100,11 +101,11 @@ c
character titleformat*130
c program version:
character*79 version
parameter(version='GREPG V2.34 plot green amplitudes')
parameter(version='GREPG V2.35 plot green amplitudes')
character*120 CVS_VERSION
parameter(CVS_VERSION=
& '$Id: grepg.f,v 1.25 2002-06-12 20:28:01 forbrig Exp $')
& '$Id: grepg.f,v 1.26 2002-09-13 21:05:02 forbrig Exp $')
c declare variables for io
character*80 filename
......@@ -158,11 +159,11 @@ c cursor routine
real pickdx
c commandline
integer maxopt, lastarg, iargc
parameter(maxopt=49)
parameter(maxopt=50)
character*3 optid(maxopt)
character*80 optarg(maxopt)
logical optset(maxopt), opthasarg(maxopt)
character*80 outdevice
character*80 outdevice,fourierfile
logical isolines, color, resolution, smooth
logical grid, scale, normalize
logical implot, replot, aliaspick, plinear, suppress, gmtoutput
......@@ -170,7 +171,7 @@ c commandline
logical deftitle, defxlab, defylab
logical phasecolor,phaseswitch,conlab,moveavg,polytrend,avgtrend
logical inccontr, setforeground, setbackground, hlsinterpol
logical omsqrscale, phasereverse
logical omsqrscale, phasereverse,fourierinput
real linelength, tracedx, charheight, contrth, contrord
real tenpower, scalingslim, normvalue, normlimit
real phaseshift,grayify, rgb_red, rgb_green, rgb_blue
......@@ -198,16 +199,17 @@ c here are the keys to our commandline options
data optid/2h-d,2h-i,2h-c,2h-r,2h-q,2h-p,2h-g,2h-l,2h-D,2h-s,2h-I,2h-R,
& 2h-P,2h-L,2h-S,2h-m,2h-G,2h-H,2h-F,2h-v,2h-N,2h-T,2h-W,2h-X,2h-Y,
& 2h-n,2h-A,2h-M,2h-C,2h-B,2h-a,2h-O,2h-f,2h-t,2h-Q,2h-K,2h-b,2h-x,
& 2h-y,2h-h,2h-e,2h-E,3h-Fc,3h-Ah,3h-Ta,3h-FC,3h-Fl,3h-Pr,3h-Ww/
& 2h-y,2h-h,2h-e,2h-E,3h-Fc,3h-Ah,3h-Ta,3h-FC,3h-Fl,3h-Pr,3h-Ww,
& 3h-ff/
data opthasarg/.TRUE.,2*.FALSE.,.TRUE.,.FALSE.,.TRUE.,.FALSE.,.TRUE.,
& 7*.FALSE.,4*.TRUE.,.FALSE.,2*.TRUE.,.FALSE.,5*.TRUE.,.FALSE.,
& 5*.TRUE.,.FALSE.,4*.TRUE.,.FALSE.,.TRUE.,.FALSE.,5*.TRUE.,
& .FALSE.,.TRUE./
& .FALSE.,2*.TRUE./
data optarg/3hx11,2*1h-,5h1.,1.,1h-,3h2.5,1h-,1h1,7*1h-,2h0.,
& 7hxyz.out,2h1.,2*1h-,5h5.,4.,5htitle,1h-,6hxlabel,6hylabel,3h-1.,
& 9hamplitude,2h0.,1h-,9hphase [],6h0.,10.,7h20,10,2,1h9,3h6,1,1h-,
& 6h4.,0.3,1h3,2*8h0.,0.,0.,1h-,2h0.,1h-,'0.,0.,0.,1','1.','NSP',
& 1h5,1h4,1h-,2h3./
& 1h5,1h4,1h-,2h3.,3hnil/
c======================================================================
c give basic information
......@@ -231,6 +233,8 @@ c give basic information
print *,
&' [-FC n] [-Fl n] [-Pr] [-Ww w]'
print *,
&' [-ff filename]'
print *,
&'or: grepg -help'
if (iargc().lt.1) stop 'ERROR: missing arguments'
call getarg(1, filename)
......@@ -389,6 +393,9 @@ c give basic information
print *,' -G file writes one-column xyz-data to file'
print *,' this file may be converted using xyz2grd.sh'
print *,' -F file read dispersion curve picks from file'
print *,' -ff filename'
print *,' read Fourier data for phasor walkout'
print *,' from file ''filename'' '
print *,' '
print *,'user interaction'
print *,'----------------'
......@@ -529,6 +536,8 @@ c campatibility to version prior V2.23
read(optarg(47), *, err=99, end=98) plpar_lscyc
phasereverse=optset(48)
read(optarg(49), *, err=99, end=98) wedgewidth
fourierinput=optset(50)
fourierfile=optarg(50)
c print *,subtitlestring
c print *,optarg(45)
......@@ -584,6 +593,14 @@ c
hints='active keys: <p,P,x,X,r,d,D,c,1,2,3,4,5,6,7,8,9,0,h,l,s,b>'
write (pickhint, 51) active_pick
c
c----------------------------------------------------------------------
c read Fourier data
if (fourierinput) then
call grepg_readfourier(fourierfile, verbose, om, nfreq, maxfreq)
else
call grepg_nofourier
endif
c----------------------------------------------------------------------
c create plot-data
c ----------------
......@@ -1494,6 +1511,7 @@ c
print *,' as created by ''gabor''.'
print *,' '
print *,' Time values (gabor matrix) are given in 1s.'
call grepg_fourierinfo
return
end
Supports Markdown
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