Commit 67e99149 authored by thomas.forbriger's avatar thomas.forbriger
Browse files

[TASK] (disper): remove obsolete 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.

disper.f is an obsolete predecessor to greda.f
parent dd7056e6
This is a legacy version of the repository. It may be incomplete as well as
inconsistent. See README.history for details. For the old stock of the
repository copyright and licence conditions apply as specified for versions
commited after 2015-03-01. Use recent versions as a base for new development.
The legacy version is only stored to keep a record of history.
# this is <Makefile>
# ----------------------------------------------------------------------------
# $Id$
#
# Copyright (c) 2008 by Thomas Forbriger (BFO Schiltach)
#
# synthetic dispersion relation
#
# REVISIONS and CHANGES
# 08/02/2008 V1.0 Thomas Forbriger
# 18/01/2011 V1.1 migrated to new scheme
#
# ============================================================================
PROGRAMS=flspherplot scanflspherfile knopoff
.PHONY: all doxydoc
all: install
.PHONY: install
install: $(addprefix $(LOCBINDIR)/,$(PROGRAMS))
$(LOCBINDIR)/%: %
mkdir -pv $(LOCBINDIR)
/bin/mv -fv $< $(LOCBINDIR)
# check mandatory environment variable settings
# ---------------------------------------------
CHECKVAR=$(if $($(1)),,$(error ERROR: missing variable $(1)))
CHECKVARS=$(foreach var,$(1),$(call CHECKVAR,$(var)))
$(call CHECKVARS,LOCINCLUDEDIR LOCLIBDIR LOCBINDIR)
FLAGS += $(MYFLAGS) -fPIC
FFLAGS += -ff2c -Wall -ffixed-line-length-0 -fno-backslash $(FLAGS)
CFLAGS += $(FLAGS)
CXXFLAGS+=-Wall $(FLAGS)
LDFLAGS+=-L$(LOCLIBDIR)
CPPFLAGS+=-I$(LOCINCLUDEDIR) $(FLAGS) -Iwosub
#----------------------------------------------------------------------
flist: Makefile $(wildcard *.inc) $(wildcard *.f) $(wildcard */*.f) \
$(wildcard *.cc) $(wildcard */*.cc */*.h)
echo $^ | tr ' ' '\n' | sort > $@
.PHONY: edit
edit: flist; vim $<
.PHONY: clean
clean:
find . -name \*.o -o -name \*.bak -o -name \*.d | \
xargs --no-run-if-empty /bin/rm -v
-/bin/rm -vf flist $(PROGRAMS)
#----------------------------------------------------------------------
# Fortran dependencies
# --------------------
%.d: %.f
echo $<: $(shell cat $< | egrep '^ +include' | cut -f 2 -d \' | sort | uniq) > $@
-include $(patsubst %.f,%.d,$(wildcard *.f))
# C++ dependencies
%.d: %.cc
$(SHELL) -ec '$(CXX) -M $(CPPFLAGS) $< \
| sed '\''s,\($(notdir $*)\)\.o[ :]*,$(dir $@)\1.o $@ : ,g'\'' \
> $@; \
[ -s $@ ] || rm -f $@'
-include $(patsubst %.cc,%.d,$(wildcard *.cc wosub/*.cc))
#----------------------------------------------------------------------
%.o: %.f
$(FC) -o $@ -c $< -O2 $(FFLAGS)
knopoff: %: %.o
$(FC) -o $@ $< -lemod -ltf $(TF_LINK_PGPLOT) -L$(LOCLIBDIR)
wosub/%.o: wosub/%.cc wosub/%.h
$(CXX) $(CXXOPT) -c -o $@ $< $(CPPFLAGS)
flspherplot: %: %.cc \
wosub/flspherclass.o wosub/flspherplotclass.o
scanflspherfile: %: %.cc
scanflspherfile flspherplot:
$(CXX) $(CXXOPT) -o $@ $^ -Wall \
-lpgplotCpp -lcpgplot -ltfxx -laff \
$(TF_LINK_PGPLOT) \
$(TF_LINK_FORTRAN) \
$(CPPFLAGS) $(LDFLAGS)
#======================================================================
# docu section
# doxygen part
# ------------
$(call CHECKVARS,TF_WWWBASEDIR TF_BROWSER)
DOXYWWWPATH=$(TF_WWWBASEDIR)/flspherplot
.PHONY: doxyclean doxyview doxydoc doxyconf
doxyclean: ;/bin/rm -rfv $(DOXYWWWPATH)
DOXYSRC=flspherplot.cc scanflspherfile.cc $(wildcard wosub/*.cc) \
$(wildcard wosub/*.h)
# create doxygen intermediate configuration
PWD=$(shell env pwd)
doxydoc.xxx: doxydoc.cfg
sed 's,<OUTPUTDIRECTORY>,$(DOXYWWWPATH),g;s,<STRIPFROMPATH>,$(PWD),g' \
$< > $@
# create commented version of doxygen configuration
doxycomm.xxx: doxydoc.cfg
/bin/cp -vf $< $@; doxygen -u $@
$(DOXYWWWPATH)/html/index.html: doxydoc.xxx $(DOXYSRC)
mkdir -vp $(DOXYWWWPATH)
doxygen $<
doxydoc: $(DOXYWWWPATH)/html/index.html
doxyview: $(DOXYWWWPATH)/html/index.html
$(TF_BROWSER) file:$< &
# ----- END OF Makefile -----
program disan
c======================================================================
c this is <disan.f>
c
c a program to do a phase-slowness-frequency
c analysis from a seismic-trace file assuming plane waves
c
c 9/1995 Thomas Forbriger
c Institut fuer Geophysik
c Universitaet Stuttgart
c
c Revisions:
c V1.0 19/09/1995 first running version
c V1.1 20/09/1995 improved with PI
c V1.2 21/09/1995 select frequency range by frequency value
c V1.3 14/02/1996 tapering of time-series
c V2.0 27/05/1997 started conversion to sff-format and fast
c processing
c
c======================================================================
c
c that's what we are going to do in subroutine
c ------------------------------ -------------
c 1. read a seismic file READDATA
c 2. stack the trace according to phase velocities DOSTACK
c 3. save stacked traces WRITESTACK
c 4. do a fourier transform FTRANS
c 5. save dispersion-array WRITEDIS
c
c======================================================================
c
c declare variables
c -----------------
c
c array dimensions
character*79 version
parameter(version='DISAN V3.0 dispersion analysis')
integer maxtraces, maxsamples, maxslo, maxslosamp, maxfreq
integer maxfree
parameter(maxtraces=100, maxsamples=100000, maxslo=80)
parameter(maxslosamp=2050, maxfreq=maxslo)
parameter(maxfree=50)
c general data
integer iargc
c seismic data
real data(maxsamples)
integer firstsample(maxtraces)
real dt(maxtraces)
real tdif(maxtraces)
real distance(maxtraces)
character*80 filename
c stack data
character*80 slodatafile
real slowness(maxslo)
real stackdata(maxslo, maxslowsamp)
integer nslo
real slodt
real maxs, mins
c fourier data
character*80 freqdatafile
complex freqdata(maxslo, maxfreq)
integer nfreq
real df, maxf
c----------------------------------------------------------------------
c
c program prolog
c
c give basic information
c ----------------------
call usage(version, maxtraces, maxslo, maxsamples,
& maxfreq, maxfree, maxslosamp)
c----------------------------------------------------------------------
c
c call subroutines
c
c read data
call readdata(maxprolog, maxtraces, maxsamples,
& maxslownesses,
& controlfile, seisdatafile,
& nprolog, ntraces, nsamples,
& seisdata, distance, dt, stmin, stsec, prolog,
& maxf, nslo, mins, maxs)
c do stacking
nslosamp=nsamples*2
c not read from controlfile:
tappercent=0.1
dotapering=.TRUE.
call dostack(maxtraces, maxsamples, maxslownesses,
& ntraces, nsamples, nslo,
& dt, stmin, stsec, mins, maxs,
& distance, seisdata, slowness, stackdata,
& dotapering, tappercent,
& maxslosamp, nslosamp)
c basename for stackfile is seisfile
write(stackdatafile(1:20),'(a20)') seisdatafile
c write information to disk
call writestack(stackdatafile,
& maxslownesses,
& nslo, mins, maxs,
& slowness, stackdata,
& maxprolog, nprolog, prolog,
& dt, version,
& maxslosamp, nslosamp)
c do the fourier transform
call ftrans(maxslownesses, maxslosamp,
& stackdata, nslo, nslosamp,
& dt, maxfreq, maxf,
& nfreq, df, freqdata)
c basename for stackfile is seisfile
write(freqdatafile(1:20),'(a20)') seisdatafile
c write information to disk
call writedis(freqdatafile,
& maxslownesses, maxfreq,
& nslo, nfreq,
& slowness, freqdata,
& maxprolog, nprolog, prolog,
& df)
stop
end
c======================================================================
c
c SUBROUTINES
c
c----------------------------------------------------------------------
c
c READDATA
c
c read seismic file
subroutine readdata(maxprolog, maxtraces, maxsamples,
& maxslownesses,
& controlfile, seisdatafile,
& nprolog, ntraces, nsamples,
& data, distance, dt, stmin, stsec, prolog,
& maxf, nslo, mins, maxs)
c declare parameters
integer maxprolog, maxtraces, maxsamples
integer maxslownesses
integer nprolog, ntraces, nsamples, nslo
real dt, stmin, stsec, mins, maxs, maxf
real distance(maxtraces), data(maxtraces, maxsamples)
character*20 controlfile, seisdatafile
character*80 prolog(maxprolog)
c local variables
c read control file
print *,'reading control file: ',controlfile
goto 1
2 stop 'ERROR: opening control file'
3 stop 'ERROR: closing control file'
4 stop 'ERROR: reading control file'
5 stop 'ERROR: unexpected end of control file'
1 open (10, file=controlfile, err=2, status='old')
read (10, *, err=4, end=5) nslo, mins, maxs, maxf
close (10, err=3)
c echo the read parameters
print *,'parameters are:'
print *,' number of slownesses: ',nslo
print *,' minimum slowness: ', mins
print *,' maximum slowness: ', maxs
print *,' up to frequency: ',maxf
c if (minf.ge.maxf) then
c print *,'ERROR: maximum frequency must be higher than ',
c & 'minimum frequency'
c stop
c endif
if (mins.ge.maxs) then
print *,'ERROR: maximum slowness must be higher than ',
& 'minimum slowness'
stop
endif
if (nslo.gt.maxslownesses) then
print *,
& 'ERROR: number of slownesses must not be higher than ',
& maxslownesses
stop
endif
c get seismic data
call seisread(seisdatafile,
& maxtraces, maxsamples, maxprolog,
& nprolog, ntraces, nsamples,
& dt, stmin, stsec,
& distance, data, prolog)
end
c----------------------------------------------------------------------
c
c DOSTACK
c
c phase velocity stack
subroutine dostack(maxtraces, maxsamples, maxslownesses,
& ntraces, nsamples, nslo,
& dt, stmin, stsec, mins, maxs,
& distance, data, slowness, stack,
& dotapering, tappercent,
& maxslosamp, nslosamp)
c declare parameters
integer maxtraces, maxsamples, maxslownesses
integer ntraces, nsamples, nslo
integer maxslosamp, nslosamp
real dt, stmin, stsec, mins, maxs
real distance(maxtraces), data(maxtraces, maxsamples)
real slowness(maxslownesses)
real stack(maxslownesses, maxslosamp)
logical dotapering
real tappercent
c declare local variables
integer trace, sample, islo, ioffset, minst, maxst
real average, ampsquare
real tapfakt, normfakt
real slostep, offset, fraction
integer midslosamp
integer winleft, winright
real tflib_costap
c tapering windows
winleft=int(float(nsamples)*tappercent)
winright=nsamples-winleft
c
c normalize data to amplitude square (energy like)
c and use a distance and time taper
c ------------------------------------------------
c
print *,'normalizing and tapering traces'
print *,
& 'trace average sqrt(amplitude**2) tapering-factor'
c go through all traces
do 1 trace=1,ntraces
c evaluate average
average=0
do 2 sample=1,nsamples
average=average+data(trace, sample)
2 continue
average=average/float(nsamples)
c remove average and evaluate amplitude square
ampsquare=0
do 3 sample=1,nsamples
data(trace, sample)=data(trace, sample)-average
ampsquare=ampsquare+data(trace, sample)**2
3 continue
c normalize to sqrt(ampsquare) and do trapezoid-tapering (if wished)
if (dotapering) then
c trapezoid-faktor
tapfakt=min((float(ntraces)*tappercent),
& float(min(trace, ntraces-trace+1)))/
& (ntraces*tappercent)
else
c no tapering
tapfakt=1
endif
c normalizing faktor
normfakt=sqrt(ampsquare)/tapfakt
c go through all samples
do 4 sample=1,nsamples
data(trace, sample)=data(trace, sample)*
& tflib_costap(sample,0,winleft,winright,nsamples)/normfakt
4 continue
c enjoy your user
write(6, 5) trace, average, sqrt(ampsquare), tapfakt
5 format(i5, f14.6, f21.7, f18.7)
1 continue
c
c stack data
c ----------
c evaluate slowness-stepwidth
slostep=(maxs-mins)/float(nslo-1)
midslosamp=int(nslosamp/2)
c go through all slownesses
do 50 islo=1,nslo
c evaluate actual slowness
slowness(islo)=mins+slostep*float(islo-1)
c do a user-feedback
print *,'stacking slowness ',
& slowness(islo),' s/m #',islo,' of ',nslo
c clear stack-array
do 52 sample=1,nsamples*2
stack(islo,sample)=0
52 continue
c go through all traces
do 51 trace=1,ntraces
c evaluate sample-shift
c traveltime=distance(trace)*slowness
c traveltime relative to beginning of trace=traveltime-(stsec+60*stmin)
c relative traveltime in samples=relative traveltime/dt
c evaluate this offset:
offset=((distance(trace)*slowness(islo))
& -(stsec+60*stmin))/dt
c offset is a real value an may lie between two samples
ioffset=int(offset)
c fraction gives us the weight of the left and the right sample to offset
fraction=offset-float(ioffset)
c evaluate the first (minst) and the last (maxst) sample to stack
c according to length of the original trace (1-nsamples) and the
c stacked trace (1-nslosamp)
minst=max(1,ioffset-midslosamp+1)
maxst=min(midslosamp-1,nslosamp+ioffset-midslosamp)
c do stacking
do 53 sample=minst,maxst
stack(islo,midslosamp-ioffset+sample)
& =stack(islo,midslosamp-ioffset+sample)
& +(float(1)-fraction)*data(trace,sample)
& +fraction*data(trace,sample+1)
53 continue
51 continue
50 continue
return
end
c----------------------------------------------------------------------
c
c WRITESTACK
c
c write stacked data
subroutine writestack(basename,
& maxslownesses,
& nslo, mins, maxs,
& slowness, stack,
& maxprolog, nprolog, prolog,
& dt, version,
& maxslosamp, nslosamp)
c declare parameters
character*80 version
character*20 basename
integer maxslownesses
integer nslo
real dt
real mins, maxs
real slowness(maxslownesses)
real stack(maxslownesses, maxslosamp)
integer maxprolog, nprolog
integer maxslosamp, nslosamp
character*80 prolog(maxprolog)
c declare local variables
integer pos
c create filename
pos=0
1 pos=pos+1
if (basename(pos:pos).ne.' ') goto 1
write(basename(pos:pos+3),'(a4)') '.slo'
print *,'writing stacked data to file ',basename
c create prolog
write(prolog(nprolog+1),'(a80)') version
write(prolog(nprolog+2),'(a8,i4,a16,f10.5,a7,f10.5,a3)')
& 'stacked ',nslo,
& ' slownesses from',
& mins,'s/m to ',maxs,'s/m'
print *,prolog(nprolog+1)
print *,prolog(nprolog+2)
nprolog=nprolog+2
c write to file
call slowwrite(basename,
& maxslownesses, maxslosamp, maxprolog,
& nprolog, nslo, nslosamp,
& dt,
& slowness, stack, prolog)
return
end
c----------------------------------------------------------------------
c
c FTRANS
c
c do a fourier transform
subroutine ftrans(maxslownesses, maxslosamp,
& stackdata, nslo, nslosamp,
& dt, maxfreq, maxf,
& nfreq, df, freqdata)
c declare parameters
integer maxslownesses, maxslosamp
real stackdata(maxslownesses, maxslosamp)
integer nslo, nslosamp
real dt, df, maxf
integer maxfreq, nfreq
complex freqdata(maxslownesses, maxfreq)
c declare local variables
complex ebase, efakt, efbase
integer islo, ifreq, sample
c evaluate frequency base-value
df=1/(float(nslosamp)*dt)
nfreq=int(maxf/df)
if (nfreq.gt.maxfreq) then
print *,
& 'ERROR: number of frequencies must not be higher than ',
& maxfreq
stop
endif
print *,'Fourier-Transform for all slownesses'
print *,nfreq,' frequencies from ',df,
& ' Hz to ',nfreq*df,' Hz'
c evaluate basevalue of exp-function
ebase=cexp((0,1)*cmplx(4*asin(float(1))/float(nslosamp)))
do 1 islo=1,nslo
efbase=(1,0)
print *,'processing slowness #',islo,' of ',nslo
do 2 ifreq=1,nfreq
efbase=efbase*ebase
efakt=(1,0)
freqdata(islo, ifreq)=(0, 0)
do 3 sample=1,nslosamp
freqdata(islo, ifreq)=freqdata(islo, ifreq)+
& efakt*cmplx(stackdata(islo, sample))
efakt=efakt*efbase
3 continue
2 continue
1 continue
return
end
c----------------------------------------------------------------------
c
c WRITEDIS
c
c write dispersion array
subroutine writedis(basename,
& maxslownesses, maxfreq,
& nslo, nfreq,
& slowness, freqdata,
& maxprolog, nprolog, prolog,
& df)
c declare parameters
character*20 basename
integer maxslownesses, maxfreq
integer nslo, nfreq
real df
real slowness(maxslownesses)
complex freqdata(maxslownesses, maxfreq)
integer maxprolog, nprolog
character*80 prolog(maxprolog)
c declare local variables
integer pos
c create filename
pos=0
1 pos=pos+1
if (basename(pos:pos).ne.' ') goto 1
write(basename(pos:pos+3),'(a4)') '.slf'
print *,'writing transformed data to file ',basename
c create prolog
write(prolog(nprolog+1),'(a22,i5,a12)')
& 'fourier-transform for ',
& nfreq,' frequencies'
write(prolog(nprolog+2),'(a4,f10.5,a7,f10.5,a3)')
& 'from', df,'Hz to ',df*nfreq,'Hz'
print *,prolog(nprolog+1)
print *,prolog(nprolog+2)
nprolog=nprolog+2
c write to file
call dispwrite(basename,
& maxslownesses, maxfreq, maxprolog,
& nprolog, nslo, nfreq,
& df,
& slowness, freqdata, prolog)
return
end
c----------------------------------------------------------------------
c
c general information
c
subroutine usage(version, maxtraces, maxslo, maxsamples,
& maxfreq, maxfree, maxslosamp)
character*(*) version
integer maxtraces, maxslo, maxsamples, maxfreq, maxfree, maxslosamp
c
character*20 param
c
print *,version
print *,'Usage: disan controlfile datafile'
print *,'or: disan -help'
call getarg(1, param)
if (param(1:5).eq.'-help') then
print *,' '
print *,' datafile name of sff-style seismic-file'
print *,' '
print *,' ns number of phase slownesses'
print *,' mins minimal phase slowness (in s/m)'
print *,' maxs maximal phase slowness (in s/m)'
print *,' maxf frequency up to which to calculate'
print *,' '
print *,' this program is compiled for:'
print *,' total number of samples: ',maxsamples
print *,' total number of traces: ',maxtraces
print *,' total number of slownesses: ',maxslo
print *,' total number of slowness samples: ',maxslosamp
print *,' total number of frequencies: ',maxfreq
print *,' total number of FREE lines: ',maxfree
stop
endif
c
return
end
c
c ----- END OF disan.f -----
program disan
c======================================================================
c this is file <disan.f>