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

added some modification options to grepg

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: 129
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 60d79bba
......@@ -12,7 +12,8 @@ F2CFLAGS=-f -u
CFLAGS=-O2
GREBOBS=grepg.o grepg_message.o grepg_dopicks.o grepg_selstyle.o \
grepg_phase.o grepg_phasewedg.o
grepg_phase.o grepg_phasewedg.o grepg_poly.o grepg_remavg.o \
grepg_contr.o
make.incdep: *.f
incdep > make.incdep
......
......@@ -55,6 +55,9 @@ c raw or smooth plots
c * introduced contour labelling
c * reorganized help section
c V2.15 16/11/99 smooth frequency dependent scaling
c V2.16 07/04/00 * added polynomial trend removal
c * and average removal
c * and contrast manipulation
c
c======================================================================
program grepg
......@@ -68,7 +71,7 @@ c
character titleformat*130
c program version:
character*79 version
parameter(version='GREPG V2.15 plot green amplitudes')
parameter(version='GREPG V2.16 plot green amplitudes')
c declare variables for io
character*80 filename
......@@ -122,7 +125,7 @@ c cursor routine
real pickdx
c commandline
integer maxopt, lastarg, iargc
parameter(maxopt=33)
parameter(maxopt=36)
character*2 optid(maxopt)
character*80 optarg(maxopt)
logical optset(maxopt), opthasarg(maxopt)
......@@ -130,11 +133,12 @@ c commandline
logical isolines, color, resolution, smooth, grid, scale, normalize
logical implot, replot, aliaspick, plinear, suppress, gmtoutput
logical readpickfile, verbose, whitepaper, deftitle, defxlab, defylab
logical phasecolor,phaseswitch,conlab,moveavg
real linelength, tracedx, charheight
logical phasecolor,phaseswitch,conlab,moveavg,polytrend,avgtrend
logical inccontr
real linelength, tracedx, charheight, contrth, contrord
real tenpower, scalingslim, normvalue, normlimit
real phaseshift,grayify
integer setlinewidth, avglength
integer setlinewidth, avglength, npolytrend
character*80 gmtxyzfile, pickfile, titlestring
character*80 xlabelstring, ylabelstring, wedgestring, phasestring
c pgplot
......@@ -154,13 +158,14 @@ c phase function pi constant
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-n,2h-A,2h-M,2h-C,2h-B,2h-a,2h-O,2h-f,2h-t,2h-Q,2h-K/
data opthasarg/.TRUE.,2*.FALSE.,.TRUE.,.FALSE.,.TRUE.,.FALSE.,.TRUE.,
& 7*.FALSE.,4*.TRUE.,.FALSE.,2*.TRUE.,.FALSE.,5*.TRUE.,.FALSE.,
& 4*.TRUE./
& 5*.TRUE.,.FALSE.,.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/
& 9hamplitude,2h0.,1h-,9hphase [],6h0.,10.,7h20,10,2,1h9,1h6,1h-,
& 6h4.,0.3/
c======================================================================
c give basic information
......@@ -176,7 +181,9 @@ c give basic information
print *,
&' [-n value] [-A label] [-M lim] [-C] [-B label]'
print *,
&' [-O i,m,s] [-a a,g] [-f n]'
&' [-O i,m,s] [-a a,g] [-f n] [-t n] [-Q}'
print *,
&' [-K o,t]'
print *,
&'or: grepg -help'
if (iargc().lt.1) stop 'ERROR: missing arguments'
......@@ -198,9 +205,14 @@ c give basic information
print *,'-------------------------'
print *,' '
print *,' -S suppress all values for zero frequency'
print *,' and values for zero slowenss'
print *,' (usefull for strange greda methods)'
print *,' -I use just imaginary part'
print *,' -R use just real part'
print *,' -t n remove polynomial trend from data using'
print *,' Legendre polynomials up to order n'
print *,' -Q remove moving average'
print *,' -K o,t increase contrast by order o and threshold t'
print *,' '
print *,' normalize data'
print *,' '
......@@ -398,6 +410,11 @@ c
moveavg=optset(33)
read(optarg(33), *) avglength
avglength=max(2,avglength)
polytrend=optset(34)
read(optarg(34), *) npolytrend
avgtrend=optset(35)
inccontr=optset(36)
read(optarg(36), *) contrord, contrth
c
c covariances
if (replot) implot=.false.
......@@ -477,6 +494,10 @@ c
do islo=1,nslo
data(islo,1)=cmplx(lowlimit)
enddo
print *,'suppress zero slowness values to lower limit ',lowlimit,'...'
do ifreq=1,nfreq
data(1,ifreq)=cmplx(lowlimit)
enddo
endif
if (optset(16)) print *,'limit scaling to slowness>=',scalingslim
c
......@@ -573,6 +594,22 @@ c
c----------------------------------------------------------------------
c
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)
endif
c
c remove average trend if requested
if (avgtrend) then
call grepg_remavg(nslo,nfreq,data,verbose, maxslo, maxfreq)
endif
c
c increase contrast
if (inccontr) then
call grepg_contr(nslo,nfreq,data,maxslo,
& maxfreq,contrth,verbose,contrord)
endif
c
c prepare maxvalue
if (plinear) then
......
c this is <grepg_contr.f>
c------------------------------------------------------------------------------
cS
c ($Source: /home/tforb/svnbuild/cvssource/CVS/thof/src/green/grepg/grepg_contr.f,v $)
c ($Id: grepg_contr.f,v 1.1 2000-04-07 16:19:54 thof Exp $)
c
c 07/04/2000 by Thomas Forbriger (IfG Stuttgart)
c
c increase contrast
c
c REVISIONS and CHANGES
c 07/04/2000 V1.0 Thomas Forbriger
c
c $Log: not supported by cvs2svn $
c
c==============================================================================
c
subroutine grepg_contr(ns,nf,d,ms,mf,t,v,o)
c
c declare parameters
real t,o
logical v
integer ns,nf,ms,mf
complex d(ms,mf)
c
cE
c declare local variables
character*(*) grepg_contr_id
parameter (grepg_contr_id='$Id: grepg_contr.f,v 1.1 2000-04-07 16:19:54 thof Exp $')
double precision maxval, x, xd
integer i,j
c
double precision f
f(o,x)=log10(1.d0+((10.d0**o)-1.d0)*x)/dble(o)
c
c------------------------------------------------------------------------------
c go
if (v) print *,'increase contrast with threshold (order: ',o,
& ' threshold: ',t,')'
do j=1,nf
maxval=0.d0
do i=1,ns
maxval=max(maxval,abs(d(i,j)))
enddo
do i=1,ns
xd=abs(d(i,j))
x=xd/maxval
if (x.gt.t) then
x=(x-t)/(1.d0-t)
x=f(o,x)*(1.d0-t)+t
else
x=-1.d0*(x-t)/t
x=t-f(o,x)*t
endif
d(i,j)=maxval*x*d(i,j)/xd
enddo
enddo
c
return
c the following line prevents the linker from removing the ID string
99 print *, grepg_contr_id
end
c
c ----- END OF grepg_contr.f -----
c this is <grepg_poly.f>
c------------------------------------------------------------------------------
cS
c ($Source: /home/tforb/svnbuild/cvssource/CVS/thof/src/green/grepg/grepg_poly.f,v $)
c ($Id: grepg_poly.f,v 1.1 2000-04-07 16:19:54 thof Exp $)
c
c 07/04/2000 by Thomas Forbriger (IfG Stuttgart)
c
c remove polynomial trend
c
c REVISIONS and CHANGES
c 07/04/2000 V1.0 Thomas Forbriger
c
c==============================================================================
c
subroutine grepg_poly(no,ns,nf,d,v, ms, mf)
c
c declare parameters
integer no, ns, nf, ms, mf
complex d(ms,mf)
logical v
c
cE
c declare local variables
character*(*) grepg_poly_id
parameter (grepg_poly_id='$Id: grepg_poly.f,v 1.1 2000-04-07 16:19:54 thof Exp $')
double precision nh, sum, msum, lp
integer i,j,n
c
c------------------------------------------------------------------------------
c go
if (v) print *, 'remove polynomial trend'
c
c here we remove a polynomial trend for marine data
n=max(0,min(6,no))
nh=dble(ns)/2.d0
do n=0,n
msum=0.d0
do j=1,nf
sum=0.d0
do i=1,ns
sum=sum+lp(n,(dble(i)/nh-1.d0))*aimag(d(i,j))
c print *,n,i,j,sum,lp(n,(dble(i)/nh-1.d0)),d(i,j),abs(d(i,j))
enddo
sum=sum/nh
do i=1,ns
d(i,j)=cmplx(real(d(i,j)),
& aimag(d(i,j))-sum*lp(n,(dble(i)/nh-1.d0)))
enddo
msum=msum+sum
enddo
do j=1,nf
sum=0.d0
do i=1,ns
sum=sum+lp(n,(dble(i)/nh-1.d0))*real(d(i,j))
c print *,n,i,j,sum,lp(n,(dble(i)/nh-1.d0)),d(i,j),abs(d(i,j))
enddo
sum=sum/nh
do i=1,ns
d(i,j)=cmplx(real(d(i,j))-sum*lp(n,(dble(i)/nh-1.d0)),
& aimag(d(i,j)))
enddo
msum=msum+sum
enddo
msum=msum/dble(nf)
if (v) print *,' polynomial order: ',n,' mean coefficient: ',msum
enddo
c
return
c the following line prevents the linker from removing the ID string
99 print *, grepg_poly_id
end
c
c----------------------------------------------------------------------
c legendre polynomials
double precision function lp(n,x)
integer n, nord
double precision x
double precision lnorm
lnorm(nord)=sqrt((2.d0*dble(nord)+1)/2.d0)
if (n.eq.0) then
lp=lnorm(0)
elseif (n.eq.1) then
lp=lnorm(1)*x
elseif (n.eq.2) then
lp=lnorm(2)*0.5d0*(3.d0*(x**2)-1.d0)
elseif (n.eq.3) then
lp=lnorm(3)*0.5d0*(5.d0*(x**3)-3.d0*x)
elseif (n.eq.4) then
lp=lnorm(4)*0.125d0*(35.d0*(x**4)-30.d0*(x**2)+3.d0)
elseif (n.eq.5) then
lp=lnorm(5)*0.125d0*(63.d0*(x**5)-70.d0*(x**3)+15.d0*x)
elseif (n.eq.6) then
lp=lnorm(6)*0.0625d0*(231.d0*(x**6)-315.d0*(x**4)+105.d0*(x**2)-5.d0)
else
stop 'illegal polynomial order!'
endif
return
end
c
c ----- END OF grepg_poly.f -----
c this is <grepg_remavg.f>
c------------------------------------------------------------------------------
cS
c ($Source: /home/tforb/svnbuild/cvssource/CVS/thof/src/green/grepg/grepg_remavg.f,v $)
c ($Id: grepg_remavg.f,v 1.1 2000-04-07 16:19:54 thof Exp $)
c
c 07/04/2000 by Thomas Forbriger (IfG Stuttgart)
c
c remove average
c
c REVISIONS and CHANGES
c 07/04/2000 V1.0 Thomas Forbriger
c
c $Log: not supported by cvs2svn $
c
c==============================================================================
c
subroutine grepg_remavg(ns,nf,d,v,ms,mf)
c
c declare parameters
integer ns,nf,ms,mf
complex d(ms,mf)
logical v
c
cE
c declare local variables
character*(*) grepg_remavg_id
parameter (grepg_remavg_id='$Id: grepg_remavg.f,v 1.1 2000-04-07 16:19:54 thof Exp $')
integer mw,i,j
parameter(mw=1000)
double precision dr(mw),di(mw),w(mw),avg,sum
c
c------------------------------------------------------------------------------
c go
c
if (mw.lt.ns) stop 'ERROR (grepg_remavg): workspace too small'
if (v) print *,'remove moving average'
avg=0.d0
do j=1,ns
do i=1,nf
dr(i)=real(d(i,j))
di(i)=aimag(d(i,j))
enddo
call rema(dr,w,nf,avg)
sum=sum+avg
call rema(di,w,nf,avg)
sum=sum+avg
do i=1,nf
d(i,j)=cmplx(dr(i),di(i))
enddo
enddo
if (v) print *,'mean average: ',sum/dble(2.d0*ns)
return
c the following line prevents the linker from removing the ID string
99 print *, grepg_remavg_id
end
c
c----------------------------------------------------------------------
c remove moving average
subroutine rema(d,w,n,sumavg)
integer n
double precision d(n), w(n), sum, ssum, pi, si, sumavg
parameter (pi=3.1415926535897931159)
integer i,j,nn
sumavg=0.d0
do i=1,n
sum=0.d0
ssum=0.d0
nn=2*i+3
do j=1,min(nn,n)
si=sin(pi*dble(j)/dble(nn))
ssum=ssum+si
sum=sum+si*d(j)
enddo
sum=sum/ssum
sumavg=sumavg+sum
w(i)=d(i)-sum
enddo
do i=1,n
d(i)=w(i)
enddo
sumavg=sumavg/dble(n)
return
end
c
c ----- END OF grepg_remavg.f -----
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