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

new handles wavenumber files

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: 967
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent ef3b150f
c this is file <grepg.f>
c======================================================================
c $Id: grepg.f,v 1.28 2002-09-17 15:36:49 forbrig Exp $
c $Id: grepg.f,v 1.29 2002-09-19 20:49:11 forbrig Exp $
c
c GREPG.F
c
......@@ -88,8 +88,9 @@ 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 V2.35a no provide phasor walkout pick
c V2.36 pick response in verbose mode
c V2.35a now provide phasor walkout pick
c V2.36 pick-response in verbose mode
c V2.37 new handle wavenumber spectra
c
c======================================================================
program grepg
......@@ -99,15 +100,15 @@ c----------------------------------------------------------------------
include 'grepg_para.inc'
include 'grepg_picks.inc'
c
logical isgrereso, isgabor
logical isgrereso, isgabor,iskspec
character titleformat*130
c program version:
character*79 version
parameter(version='GREPG V2.36 plot green amplitudes')
parameter(version='GREPG V2.37 plot green amplitudes')
character*120 CVS_VERSION
parameter(CVS_VERSION=
& '$Id: grepg.f,v 1.28 2002-09-17 15:36:49 forbrig Exp $')
& '$Id: grepg.f,v 1.29 2002-09-19 20:49:11 forbrig Exp $')
c declare variables for io
character*80 filename
......@@ -555,12 +556,22 @@ c
c get datafile
call greenread(filename, debug,
& maxslo, maxfreq, slo, om,
& nslo, nfreq, data, weight, isgrereso, isgabor)
& nslo, nfreq, data, weight, isgrereso, isgabor, iskspec)
c check dimensions for quick-mode
print *,' '
print *,'dimensions of dataset read:'
print *,' slowness: ',nslo
print *,' frequency: ',nfreq
if (iskspec) then
print *,' wavenumber: ',nslo
elseif (isgabor) then
print *,' time: ',nslo
else
print *,' slowness: ',nslo
endif
if (isgrereso) then
print *,' t-slowness: ',nfreq
else
print *,' frequency: ',nfreq
endif
c
pi2=8.*atan(1.)
c report
......@@ -582,8 +593,12 @@ c report
maxs=mins+du*(nslo-1)
if (isgabor) then
print *,' time range: ',mins,'s - ',maxs,'s dt=',du,'s'
elseif (iskspec) then
print *,' wavenumber range: ',mins,'1/km - ',maxs,
& '1/km dk=',du,'1/km'
else
print *,' slowness range: ',mins,'s/km - ',maxs,'s/km du=',du,'s/km'
print *,' slowness range: ',mins,'s/km - ',maxs,
& 's/km du=',du,'s/km'
endif
print *,' '
c----------------------------------------------------------------------
......@@ -628,20 +643,25 @@ c
enddo
endif
if (suppress) then
print *,'suppress zero frequency values to lower limit ',lowlimit,'...'
print *,'suppress zero frequency values to lower limit ',
& lowlimit,'...'
do islo=1,nslo
data(islo,1)=cmplx(lowlimit)
enddo
print *,'suppress zero slowness values to lower limit ',lowlimit,'...'
print *,'suppress zero slowness values to lower limit ',
& lowlimit,'...'
do ifreq=1,nfreq
data(1,ifreq)=cmplx(lowlimit)
enddo
endif
print *,' limit rescaling factor evaluation to data with slowness>=',
print *,
& ' limit rescaling factor evaluation to data with slowness>=',
& scalingslim
print *,' limit maximum amplitude plot scale to data with slowness>='
print *,
& ' limit maximum amplitude plot scale to data with slowness>='
& ,maxilim
print *,'limit evaluation of normalizing factor to data with slowness>=',
print *,
& 'limit evaluation of normalizing factor to data with slowness>=',
& normlimit
c
c scale each frequency to the same maximum amplitude
......@@ -691,7 +711,8 @@ c endif
maxvalue=max(abs(data(nslo,ifreq)),lowlimit)
do islo=1,nslo
slp=mins+du*(islo-1)
if (slp.ge.scalingslim) maxvalue=max(maxvalue,abs(data(islo,ifreq)))
if (slp.ge.scalingslim) maxvalue=
& max(maxvalue,abs(data(islo,ifreq)))
enddo
scaling(ifreq)=1./maxvalue
enddo
......@@ -752,7 +773,8 @@ c
do ifreq=1,nfreq
do islo=1,nslo
slp=mins+du*(islo-1)
if (slp.ge.normlimit) maxvalue=max(maxvalue,abs(data(islo,ifreq)))
if (slp.ge.normlimit) maxvalue=
& max(maxvalue,abs(data(islo,ifreq)))
enddo
enddo
maxvalue=maxvalue/max(normvalue,lowlimit)
......@@ -769,8 +791,8 @@ c here we start data modification and extraction
c
c remove polynomial trend if requested
if (polytrend) then
call grepg_poly(npolytrend,nslo,nfreq,data,verbose, maxslo, maxfreq,
& npolymode)
call grepg_poly(npolytrend,nslo,nfreq,data,verbose,
& maxslo, maxfreq, npolymode)
endif
c
c remove average trend if requested
......@@ -860,7 +882,8 @@ c come with logarithmic scales anyway.
print *,frp,slp,val1,val2,val3,val4,value
print *,max(lowlimit,abs(data(rslo,rfre))),rslo,rfre
print *,max(lowlimit,abs(data(rslo,rfre+1))),rslo,rfre+1
print *,max(lowlimit,abs(data(rslo+1,rfre+1))),rslo+1,rfre+1
print *,max(lowlimit,abs(data(rslo+1,rfre+1))),
& rslo+1,rfre+1
print *,max(lowlimit,abs(data(rslo+1,rfre))),rslo+1,rfre
endif
c phase values
......@@ -917,7 +940,8 @@ c amplitude
plotdata(islo,ifreq)=value
c phase
plotphase(islo,ifreq)=
& atan2(aimag(data(islo,ifreq)),real(data(islo,ifreq)))*180./pi
& atan2(aimag(data(islo,ifreq)),
& real(data(islo,ifreq)))*180./pi
if (plotphase(islo,ifreq).lt.0.)
& plotphase(islo,ifreq)=plotphase(islo,ifreq)+360.
if (phasereverse) plotphase(islo,ifreq)=
......@@ -992,6 +1016,10 @@ c
titleformat='("gabor-matrix (spectrogram) - plot of ",a,a)'
xlabel='frequency (Hz)'
ylabel='time (s)'
elseif (iskspec) then
titleformat='("wavenumber/frequency - plot of ",a,a)'
xlabel='frequency (Hz)'
ylabel='wavenumber (1/km)'
else
titleformat='("phase-slowness/frequency - plot of ",a,a)'
xlabel='frequency (Hz)'
......@@ -1350,7 +1378,10 @@ c 01/07/99 allow taper files
c
subroutine greenread(filename, debug,
& maxslo, maxfreq, slo, om,
& nslo, nom, green, weight, isgrereso, isgabor)
& nslo, nom, green, weight, isgrereso, isgabor,
& iskspec)
c
include 'grepg_dim.inc'
c
character filename*(*)
logical debug
......@@ -1361,18 +1392,14 @@ c complex green(maxfreq, maxslo)
real om(maxfreq), slo(maxslo)
integer nslo, nom
integer inmagic, cpu, match
character*4 cmagic, incmagic, wcmagic, pcmagic, gcmagic
parameter(cmagic='1234')
parameter(wcmagic='123S')
parameter(pcmagic='123P')
parameter(gcmagic='123G')
equivalence(incmagic, inmagic)
integer lu
parameter(lu=20)
integer i,j
logical isgreen,isgrereso,isgabor
logical isgreen,isgrereso,isgabor,iskspec
character*4 incmagic
real pi2
parameter(pi2=2.*3.1415926535897931159979634685)
parameter(pi2=2.*3.141592653589)
c
print *,'read green file ',filename(1:index(filename, ' '))
open(lu, file=filename, form='unformatted', status='old', err=99)
......@@ -1380,6 +1407,7 @@ c
c check byte sex
isgrereso=.false.
isgabor=.false.
iskspec=.false.
read(lu, err=98, end=97) inmagic
call tf_bytesex(cmagic, inmagic, cpu, match)
if (cpu.eq.1) then
......@@ -1403,6 +1431,7 @@ c check byte sex
print *,'matching bytesex (grereso output) - good...'
isgreen=.true.
isgrereso=.true.
iskspec=.false.
elseif (match.eq.2) then
print *,'bytesex not matching (grereso output)',
& ' - we will have to swap!'
......@@ -1410,6 +1439,7 @@ c check byte sex
else
isgreen=.false.
isgrereso=.false.
iskspec=.false.
print *,'bytesex read (',incmagic,') ist not grereso green'
call tf_bytesex(gcmagic, inmagic, cpu, match)
if (match.eq.1) then
......@@ -1421,9 +1451,27 @@ c check byte sex
& ' - we will have to swap!'
stop 'ERROR: do not know how to do that...'
else
print *,'bytesex read (',incmagic,') ist not gabor green'
isgreen=.false.
isgabor=.false.
isgrereso=.false.
iskspec=.false.
print *,'bytesex read (',incmagic,') ist not grereso green'
call tf_bytesex(kcmagic, inmagic, cpu, match)
if (match.eq.1) then
print *,'matching bytesex (wavenumber spectrum) - good...'
isgreen=.true.
isgabor=.false.
iskspec=.true.
elseif (match.eq.2) then
print *,'bytesex not matching (wavenumber spectrum)',
& ' - we will have to swap!'
stop 'ERROR: do not know how to do that...'
else
print *,'bytesex read (',incmagic,
& ') ist not wavenumber green'
isgreen=.false.
isgabor=.false.
iskspec=.false.
endif
endif
endif
endif
......@@ -1494,11 +1542,9 @@ c
c comments on the file format
c
subroutine grepg_filecomment
character*4 cmagic, wcmagic, pcmagic, gcmagic
parameter(cmagic='1234')
parameter(wcmagic='123S')
parameter(pcmagic='123P')
parameter(gcmagic='123G')
c
include 'grepg_dim.inc'
c
print *,' '
print *,'File Formats used by grepg:'
print *,' '
......@@ -1522,10 +1568,21 @@ c
print *,'p,p-resolution (magic number: ',pcmagic,')'
print *,' as created by ''grereso''.'
print *,' '
print *,' In the case of grereso files verbose messages'
print *,' regarding frequency concern the test-slowness value.'
print *,' '
print *,'Gabor matrix (magic number: ',gcmagic,')'
print *,' as created by ''gabor''.'
print *,' '
print *,' In the case of gabor files verbose messages regarding'
print *,' phase slowness concern time.'
print *,' Time values (gabor matrix) are given in 1s.'
print *,' '
print *,'Wavenumber spectrum (magic number: ',kcmagic,')'
print *,' as used by ''flgevas''.'
print *,' '
print *,' In the case of wavenumber Green files verbose'
print *,' messages regarding slowness concern wavenumber.'
call grepg_fourierinfo
return
end
......
c this is <grepg_dim.inc>
c------------------------------------------------------------------------------
c $Id: grepg_dim.inc,v 1.5 2002-05-03 18:12:33 forbrig Exp $
c $Id: grepg_dim.inc,v 1.6 2002-09-19 20:49:12 forbrig Exp $
c
c 25/11/98 by Thomas Forbriger (IfG Stuttgart)
c
......@@ -10,6 +10,7 @@ c REVISIONS and CHANGES
c 25/11/98 V1.0 Thomas Forbriger
c 21/01/02 V1.1 now use different dimension for read picks
c is useful for flspherplot
c 19/09/02 V1.2 provide magic numbers - for coherency
c
c==============================================================================
c
......@@ -17,4 +18,22 @@ c
parameter(maxpicktypes=10, maxpicks=400, maxpickdim=7*maxpicktypes)
c
c magic numbers
character*4 cmagic, wcmagic, pcmagic, gcmagic, kcmagic
c
c original green spectra
parameter(cmagic='1234')
c
c data taper weights for gremlin (real numbers)
parameter(wcmagic='123S')
c
c p,p-resolution as created by grereso
parameter(pcmagic='123P')
c
c gabor matrix as created by gabor
parameter(gcmagic='123G')
c
c wavenumber spectrum as taken from flgevas
parameter(kcmagic='123K')
c
c ----- END OF grepg_dim.inc -----
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