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

nice program; brand new

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: 2280
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent aea73ab4
#
# $Id: Makefile,v 1.7 2006-06-02 09:40:26 tforb Exp $
# $Id: Makefile,v 1.8 2007-05-11 09:26:01 tforb Exp $
#
# Makefile for src/green/disan
#
......@@ -44,7 +44,7 @@ greda: greda.o greda_phasor.o
-L$(LOCLIBDIR) -L$(SERVERLIBDIR)
newprog greda
hivgabor gabor phadi: %: %.o
hivexpanco hivgabor gabor phadi: %: %.o
gcc -o $@ $< -lsffu -lts -ltf -llapex -lblasex \
-lnumrec -ltime -lsff $(PGPLOTLIB) $(F2CLIB) \
-L$(LOCLIBDIR) -L$(SERVERLIBDIR)
......
c this is <hivexpanco.f>
c ----------------------------------------------------------------------------
c ($Id: hivexpanco.f,v 1.1 2007-05-11 08:58:29 tforb Exp $)
c ($Id: hivexpanco.f,v 1.2 2007-05-11 09:26:01 tforb Exp $)
c
c Copyright (c) 2007 by Thomas Forbriger (BFO Schiltach)
c
......@@ -12,80 +12,14 @@ c
c ============================================================================
c
program hivexpanco
c
character*(*) version
parameter(version='HIVEXPANCO V1.0 calculate H/V for Fourier Bessel expansion coefficients')
character*(*) HIVEXPANCO_CVS_ID
parameter(HIVEXPANCO_CVS_ID='$Id: hivexpanco.f,v 1.1 2007-05-11 08:58:29 tforb Exp $')
c
c commandline
integer maxopt, lastarg, iargc
character*80 argument
parameter(maxopt=2)
character*2 optid(maxopt)
character*40 optarg(maxopt)
logical optset(maxopt), opthasarg(maxopt)
c debugging
logical debug, verbose
c here are the keys to our commandline options
data optid/2h-d, 2h-v/
data opthasarg/2*.FALSE./
data optarg/2*1h-/
c
c------------------------------------------------------------------------------
c basic information
c
c
argument=' '
if (iargc().eq.1) call getarg(1, argument)
if ((argument(1:5).eq.'-help').or.(iargc().lt.1)) then
print *,version
print *,'Usage: hivexpanco arguments'
print *,' or: hivexpanco -help'
if (argument(1:5).ne.'-help')
& stop 'ERROR: wrong number of arguments'
print *,' '
print *,'j'
print *,' '
print *,HIVEXPANCO_CVS_ID
stop
endif
c
c------------------------------------------------------------------------------
c read command line arguments
c
call tf_cmdline(1, lastarg, maxopt, optid,
& optarg, optset, opthasarg)
debug=optset(1)
verbose=optset(2)
c
c------------------------------------------------------------------------------
c go
c
stop
end
c
c this is <hivgabor.f>
c ----------------------------------------------------------------------------
c ($Id: hivexpanco.f,v 1.1 2007-05-11 08:58:29 tforb Exp $)
c
c Copyright (c) 2006 by Thomas Forbriger (BFO Schiltach)
c
c Calculate H/V from gabor results
c
c REVISIONS and CHANGES
c 02/06/2006 V1.0 Thomas Forbriger
c
c ============================================================================
c
program hivgabor
c
character*(*) version
parameter(version=
& 'HIVGABOR V1.0 Calculate H/V from gabor results')
character*(*) HIVGABOR_CVS_ID
parameter(HIVGABOR_CVS_ID=
& '$Id: hivexpanco.f,v 1.1 2007-05-11 08:58:29 tforb Exp $')
& 'HIVEXPANCO V1.0'//
& ' calculate H/V for Fourier Bessel expansion coefficients')
character*(*) HIVEXPANCO_CVS_ID
parameter(HIVEXPANCO_CVS_ID=
& '$Id: hivexpanco.f,v 1.2 2007-05-11 09:26:01 tforb Exp $')
c
c datasets
character*80 hinfile, vinfile, outfile
......@@ -94,21 +28,21 @@ c datasets
c magic number for binary file identification
integer magic
character*4 cmagic
parameter(cmagic='123G')
parameter(cmagic='1234')
c greens function
integer maxsamples,maxom
parameter(maxsamples=3000,maxom=maxsamples)
complex hgabor(maxsamples, maxom)
complex vgabor(maxsamples, maxom)
complex hvgabor(maxsamples, maxom)
real htime(maxsamples), hom(maxom)
real vtime(maxsamples), vom(maxom)
real hivtime(maxsamples), hivom(maxom)
integer maxp,maxom
parameter(maxp=3000,maxom=maxp)
complex hgreen(maxp, maxom)
complex vgreen(maxp, maxom)
complex hvgreen(maxp, maxom)
real htime(maxp), hom(maxom)
real vtime(maxp), vom(maxom)
real hivtime(maxp), hivom(maxom)
c
c processing
integer hnsamples,hnom
integer vnsamples,vnom
integer hivnsamples,hivnom
integer hnp,hnom
integer vnp,vnom
integer hivnp,hivnom
logical overwrite, waterlevel
real water, rms
c
......@@ -135,9 +69,9 @@ c
if (iargc().eq.1) call getarg(1, argument)
if ((argument(1:5).eq.'-help').or.(iargc().lt.3)) then
print *,version
print *,'Usage: hivgabor [-v] [-o] [-w val]'
print *,' Hinfile Vinfile outfile'
print *,' or: hivgabor -help'
print *,'Usage: hivexpanco [-v] [-o] [-w val]'
print *,' Hinfile Vinfile outfile'
print *,' or: hivexpanco -help'
if (argument(1:5).ne.'-help')
& stop 'ERROR: wrong number of arguments'
print *,' '
......@@ -145,17 +79,16 @@ c
print *,'Vinfile vertical component input file'
print *,'outfile output file name'
print *,' '
print *,'The program calculates the complex ratio of gabor'
print *,'H/V from a horizontal component and vertical component'
print *,'gabor matrix.'
print *,'All files are in the format produced by gabor.f'
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 *,'-v be verbose'
print *,'-o overwrite existing output file'
print *,'-w val water level for denominator as a fraction'
print *,' of signal rms for a given time'
print *,' '
print *,HIVGABOR_CVS_ID
print *,HIVEXPANCO_CVS_ID
stop
endif
c
......@@ -176,15 +109,15 @@ c
c
c------------------------------------------------------------------------------
c go
call gaborread(hinfile, verbose, maxsamples, maxom,
& htime, hom, hnsamples, hnom, hgabor)
call gaborread(vinfile, verbose, maxsamples, maxom,
& vtime, vom, vnsamples, vnom, vgabor)
call greenread(hinfile, verbose, maxp, maxom,
& htime, hom, hnp, hnom, hgreen)
call greenread(vinfile, verbose, maxp, maxom,
& vtime, vom, vnp, vnom, vgreen)
c
c check consistency
if (vnsamples.ne.hnsamples)
if (vnp.ne.hnp)
& stop 'ERROR: inconsistent number of samples!'
hivnsamples=hnsamples
hivnp=hnp
if (vnom.ne.hnom)
& stop 'ERROR: inconsistent number of samples!'
hivnom=hnom
......@@ -198,7 +131,7 @@ c check consistency
endif
hivom(i)=hom(i)
enddo
do i=1,hnsamples
do i=1,hnp
if (htime(i).lt.1.e-20) then
if (abs(htime(i)-vtime(i)).gt.1.e-20)
& stop 'ERROR: inconsistent frequencies!'
......@@ -212,47 +145,47 @@ c
c calculate ratio
if (verbose.and.waterlevel)
& print *,'use waterlevel of ',water,' *rms'
do j=1,hivnsamples
do j=1,hivnp
if (waterlevel) then
rms=0.
do i=1,hivnom
rms=rms+abs(hgabor(j,i)*conjg(hgabor(j,i)))
rms=rms+abs(vgabor(j,i)*conjg(vgabor(j,i)))
rms=rms+abs(hgreen(j,i)*conjg(hgreen(j,i)))
rms=rms+abs(vgreen(j,i)*conjg(vgreen(j,i)))
enddo
rms=sqrt(rms)/(2.*hivnom)
do i=1,hivnom
hvgabor(j,i)=hgabor(j,i)/(vgabor(j,i)+water*rms)
hvgreen(j,i)=hgreen(j,i)/(vgreen(j,i)+water*rms)
enddo
else
do i=1,hivnom
hvgabor(j,i)=hgabor(j,i)/vgabor(j,i)
hvgreen(j,i)=hgreen(j,i)/vgreen(j,i)
enddo
endif
enddo
c
c write gabor file (easy to use)
c write green file (easy to use)
c
if (overwrite) then
if (verbose) print *,'opening gabor file ',
if (verbose) print *,'opening green file ',
& outfile(1:index(outfile,' ')),' - overwrite mode'
open(lu, file=outfile, form='unformatted', err=98)
else
if (verbose) print *,'opening gabor file ',
if (verbose) print *,'opening green file ',
& outfile(1:index(outfile,' '))
open(lu, file=outfile, status='new', form='unformatted', err=98)
endif
call tf_magic(cmagic, magic)
write(lu, err=97) magic
write(lu, err=97) hivnom, hivnsamples
write(lu, err=97) hivnom, hivnp
write(lu, err=97) (hivom(i), i=1,hivnom),
& (hivtime(i), i=1,hivnsamples)
write(lu, err=97) ((hvgabor(j,i), i=1,hivnom), j=1,hivnsamples)
& (hivtime(i), i=1,hivnp)
write(lu, err=97) ((hvgreen(j,i), i=1,hivnom), j=1,hivnp)
close(lu, err=96)
c
stop
98 stop 'ERROR: opening gabor file'
97 stop 'ERROR: writing gabor file'
96 stop 'ERROR: closing gabor file'
98 stop 'ERROR: opening green file'
97 stop 'ERROR: writing green file'
96 stop 'ERROR: closing green file'
end
c
c======================================================================
......@@ -261,23 +194,23 @@ c this routine reads in the data
c
c 01/07/99 allow taper files
c
subroutine gaborread(filename, verbose,
& maxsamples, maxfreq, time, om,
& nsamples, nom, gabor)
subroutine greenread(filename, verbose,
& maxp, maxfreq, slo, om,
& np, nom, green)
c
character filename*(*)
logical verbose
integer maxsamples, maxfreq
complex gabor(maxsamples, maxfreq)
real om(maxfreq), time(maxsamples)
integer nsamples, nom
integer maxp, maxfreq
complex green(maxp, maxfreq)
real om(maxfreq), slo(maxp)
integer np, nom
integer inmagic, cpu, match
equivalence(incmagic, inmagic)
integer lu
parameter(lu=20)
integer i,j
character*4 incmagic, cmagic
parameter(cmagic='123G')
parameter(cmagic='1234')
real pi2
parameter(pi2=2.*3.141592653589)
c
......@@ -302,26 +235,24 @@ c check byte sex
stop 'ERROR: do not know how to do that...'
else
print *,'magic number read (',incmagic,
& ') does not indicate gabor file!'
& ') does not indicate green file!'
stop 'ERROR: wrong file type!'
endif
c
read(lu, err=98, end=97) nom, nsamples
if ((nom.gt.maxfreq).or.(nsamples.gt.maxsamples)) then
read(lu, err=98, end=97) nom, np
if ((nom.gt.maxfreq).or.(np.gt.maxp)) then
close(lu, err=96)
stop 'ERROR: data exceeds array bounds'
endif
read(lu, err=98, end=97) (om(i), i=1,nom), (time(i), i=1,nsamples)
read(lu, err=98, end=97) ((gabor(j,i), i=1,nom), j=1,nsamples)
read(lu, err=98, end=97) (om(i), i=1,nom), (slo(i), i=1,np)
read(lu, err=98, end=97) ((green(j,i), i=1,nom), j=1,np)
close(lu, err=96)
print *,'gabor file read and closed'
print *,'green file read and closed'
return
99 stop 'ERROR: opening gabor file'
98 stop 'ERROR: reading gabor file'
97 stop 'ERROR: reading gabor file - unexpected end'
96 stop 'ERROR: closing gabor file'
99 stop 'ERROR: opening green file'
98 stop 'ERROR: reading green file'
97 stop 'ERROR: reading green file - unexpected end'
96 stop 'ERROR: closing green file'
end
c ----- END OF hivgabor.f -----
c
c ----- END OF hivexpanco.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