Commit 67b9c79d authored by Simone Butzer's avatar Simone Butzer

initial commit

source code and folder structure
parents
# this is <Makefile>
# ----------------------------------------------------------------------------
# ($Id: Makefile 4001 2011-06-10 12:10:09Z lrehor $)
#
# 25/10/2000 by Thomas Forbriger (IMGF Frankfurt)
#
# libseife Makefile
#
# This library contains code which was extracted from the program seife.f
# by Erhard Wielandt. The original version was written in 1984 at ETH Zurich.
# The code can be obtained through
# http://www.software-for-seismometry.de/
#
# The code from seife.f was extracted to a library by Wolfgang Friederich in
# 1996. The current version of libseife is a an extended code collection based
# on the abovementioned code.
#
# REVISIONS and CHANGES
# 25/10/2000 V1.0 Thomas Forbriger
# 14/12/2007 V1.1 g77 compilation is the default now
# 17/12/2007 V1.2 moved to gfortran
# 18/09/2010 V1.3 start migration to new SVN scheme
# discard f2c option (fallback)
# check variables
# 15/11/2010 V1.4 do not use fdep.sh
# 17/11/2010 V1.5 libseifemon.a is out of use
# 17/01/2011 V1.6 distinguish library creation and installation
# 10/06/2011 V1.7 preparing Makefile for export of cseife to
# the DENISE code by Lisa Rehor
#
# ============================================================================
#
.PHONY: all
all: install doc
.PHONY: doc
doc:
LIBRARIES=libcseife.a
.PHONY: install
install:
$(MAKE) $(LIBRARIES)
# ============================================================================
CLIBSRC=$(wildcard cseife*.c)
HEADERS=$(shell find . -name \*.h)
.PHONY: clean edit
clean:
-find . -name \*.d | xargs --no-run-if-empty /bin/rm -v
-find . -name \*.o | xargs --no-run-if-empty /bin/rm -v
-find . -name \*.bak | xargs --no-run-if-empty /bin/rm -v
-/bin/rm -fv flist
-/bin/rm -fv *.a
edit: flist ; vim flist
flist: $(wildcard Makefile *.c *.h)
echo $^ | tr ' ' '\n' | sort > $@
AR=ar
AS=as
RANLIB=ranlib
CC=mpiicc
FLAGS += $(MYFLAGS)
CFLAGS += -O3 -w $(FLAGS)
CPPFLAGS=$(addprefix -I,.) \
$(FLAGS)
.c.o:
$(CC) $(CFLAGS) $< -c $(CPPFLAGS)
#======================================================================
# reinstall target
# is necessary in case of header file problems (e.g. remcmmnt not installed)
.PHONY: reinstall
reinstall:
$(MAKE) clean
$(MAKE) install
#======================================================================
# pattern rules
# -------------
%.d: %.c
$(SHELL) -ec '$(CC) -M $(CPPFLAGS) $< \
| sed '\''s,\($(notdir $*)\)\.o[ :]*,$(dir $@)\1.o $@ : ,g'\'' \
> $@; \
[ -s $@ ] || rm -f $@'
#======================================================================
# library part
# ------------
include $(patsubst %.c,%.d,$(CLIBSRC))
libcseife.a: $(patsubst %.c,%.o,$(CLIBSRC))
libcseife.a:
$(AR) rcv $@ $^
$(RANLIB) $@
#
# ----- END OF Makefile -----
this is <README>
============================================================================
cseife
------
$Id$
============================================================================
This directory contains source code for the library cseife. It was exportet
from TF_Software to DENISE and 3DelasticFWI. DENISE is hosted in the Trac-/svn-project
FWI_elastic (http://gpitrsvn.gpi.uni-karlsruhe.de/repos/FWI_elastic/DENISE).
3DFWI is hosted in the Trac-/svn-project
FWI_elastic (http://gpitrsvn.gpi.uni-karlsruhe.de/repos/FWI_elastic/3DelasticFWI).
The copyright of the source code is held by different persons:
* cseife.c, cseife.h:
Copyright (c) 2005 by Thomas Forbriger (BFO Schiltach)
* cseife_deriv.c, cseife_gauss.c, cseife_rekfl.c, cseife_rfk.c and
cseife_tides.c:
Copyright 1984 by Erhard Wielandt
This algorithm was part of seife.f. A current version of seife.f can be
obtained from http://www.software-for-seismometry.de/
Exports:
Each time this source code is exported to DENISE and 3DFWI the current version is
tagged. The directory containing the tagged version is named by the day
of the export (yyyymmdd) and can be found in
http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/export/cseife.
10/06/2011 First export to Denise of the version tagged in 20110610
20/01/2012 First export to 3DelasticFWI of the version tagged in 20120120
----- END OF README -----
/*! \file cseife.c
* \brief make seife functions available within C code (implementation)
*
* ----------------------------------------------------------------------------
*
* $Id: cseife.c 3381 2010-11-15 12:33:00Z tforb $
* \author Thomas Forbriger
* \date 14/01/2005
*
* make seife functions available within C code (implementation)
*
* Copyright (c) 2005 by Thomas Forbriger (BFO Schiltach)
*
* REVISIONS and CHANGES
* - 14/01/2005 V1.0 Thomas Forbriger
* - 11/07/2005 V1.1 support debug mode
* - 15/11/2010 V1.2 do not use tfmacros.h
*
* ============================================================================
*/
#define TF_CSEIFE_C_VERSION \
"TF_CSEIFE_C V1.2"
#define TF_CSEIFE_C_CVSID \
"$Id: cseife.c 3381 2010-11-15 12:33:00Z tforb $"
#include <cseife.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* initialize global flags */
struct seife_flags seife_global_flags={ 0 };
/* enter debug mode */
void seife_debug_mode_on()
{
seife_global_flags.debug=1;
}
/* Butterworth lowpass (period t0, order o) */
void seife_lpb(double* x, int n, double dt, double t0, int o)
{
seife_butterworth(x, n, dt, IIRlpb, t0, o);
}
/* Butterworth highpass (period t0, order o) */
void seife_hpb(double* x, int n, double dt, double t0, int o)
{
seife_butterworth(x, n, dt, IIRhpb, t0, o);
}
/* 2nd order lowpass (period t0, damping h) */
void seife_lp2(double* x, int n, double dt, double t0, double h)
{
seife_rekfl(x, x, n, seife_rfk(IIRlp2, t0/dt, h, 0., 0.));
}
/* 2nd order highpass (period t0, damping h) */
void seife_hp2(double* x, int n, double dt, double t0, double h)
{
seife_first(x, n);
seife_rekfl(x, x, n, seife_rfk(IIRhp2, t0/dt, h, 0., 0.));
}
/* 2nd order bandpass (period t0, damping h) */
void seife_bp2(double* x, int n, double dt, double t0, double h)
{
seife_rekfl(x, x, n, seife_rfk(IIRbp2, t0/dt, h, 0., 0.));
}
/* 1st order lowpass (period t0) */
void seife_lp1(double* x, int n, double dt, double t0)
{
seife_rekfl(x, x, n, seife_rfk(IIRlp1, t0/dt, 0., 0., 0.));
}
/* 1st order highpass (period t0) */
void seife_hp1(double* x, int n, double dt, double t0)
{
seife_first(x, n);
seife_rekfl(x, x, n, seife_rfk(IIRhp1, t0/dt, 0., 0., 0.));
}
/* integration (time constant t0) */
void seife_int(double* x, int n, double dt, double t0)
{
double t00=t0;
if (t00 < 1.e-22) { t00=1.; }
seife_rekfl(x, x, n, seife_rfk(IIRint, t00/dt, 0., 0., 0.));
}
/* 1st order highpass equalizer (former period t0s, new period t0) */
void seife_he1(double* x, int n, double dt, double t0s, double t0)
{
seife_rekfl(x, x, n, seife_rfk(IIRhe1, t0/dt, 0., t0s/dt, 0.));
}
/* 1st order lowpass equalizer (former period t0s, new period t0) */
void seife_le1(double* x, int n, double dt, double t0s, double t0)
{
struct seife_IIRcoef coef=seife_rfk(IIRle1, t0/dt, 0., t0s/dt, 0.);
double fac=t0s/t0;
coef.f0 *= fac;
coef.f1 *= fac;
seife_rekfl(x, x, n, coef);
}
/* 2nd order highpass equalizer (former: period t0s and damping hs,
* new: period t0 and damping h)
*/
void seife_he2(double* x, int n, double dt,
double t0s, double hs, double t0, double h)
{
if (seife_global_flags.debug)
{
fprintf(stderr,
"DEBUG (seife_he2): n=%d, dt=%f, t0s=%f, hs=%f, t0=%f, h=%f\n",
n, dt, t0s, hs, t0, h);
}
seife_rekfl(x, x, n, seife_rfk(IIRhe2, t0/dt, h, t0s/dt, hs));
}
/* 2nd order lowpass equalizer (former: period t0s and damping hs,
* new: period t0 and damping h)
*/
void seife_le2(double* x, int n, double dt,
double t0s, double hs, double t0, double h)
{
struct seife_IIRcoef coef=seife_rfk(IIRle2, t0/dt, h, t0s/dt, hs);
double fac=t0s/t0;
fac *= fac;
coef.f0 *= fac;
coef.f1 *= fac;
coef.f2 *= fac;
seife_rekfl(x, x, n, coef);
}
/* set baseline to first value */
void seife_first(double* x, int n)
{
double x0=x[0];
int i;
for (i=0; i<n; ++i) { x[i] -= x0; }
}
/* general Butterworth filter framework */
void seife_butterworth(double* x, int n, double dt,
enum seife_filter_type IIRtype,
double t0, int m)
{
enum seife_filter_type t1st, t2nd;
if (IIRtype == IIRhpb)
{
seife_first(x, n);
t1st=IIRhp1;
t2nd=IIRhp2;
}
else if (IIRtype == IIRlpb)
{
t1st=IIRlp1;
t2nd=IIRlp2;
}
else
{
SEIFE_CHECKERROR( 1, "seife_butterworth", "illegal IIRtype" )
}
int mm=m/2;
if (m > (2*mm))
{ seife_rekfl(x, x, n, seife_rfk(t1st, t0/dt, 0., 0., 0.)); }
double pih=2.*atan(1.);
int j;
for (j=1; j<=mm; ++j)
{
double h=sin(pih*(2*j-1)/m);
seife_rekfl(x, x, n, seife_rfk(t2nd, t0/dt, h, 0., 0.));
}
}
/* ----- END OF cseife.c ----- */
/*! \file cseife.h
* \brief make seife functions available within C code (prototypes)
*
* ----------------------------------------------------------------------------
*
* $Id: cseife.h 3381 2010-11-15 12:33:00Z tforb $
* \author Thomas Forbriger
* \date 14/01/2005
*
* make seife functions available within C code (prototypes)
*
* Copyright (c) 2005 by Thomas Forbriger (BFO Schiltach)
*
* REVISIONS and CHANGES
* - 14/01/2005 V1.0 Thomas Forbriger
* - 11/07/2005 V1.1 introduce global flags for debug mode
* - 15/11/2010 V1.2 provide CHECKERROR macro
*
* ============================================================================
*/
/* include guard */
#ifndef TF_CSEIFE_H_VERSION
#define TF_CSEIFE_H_VERSION \
"TF_CSEIFE_H V1.2"
#define TF_CSEIFE_H_CVSID \
"$Id: cseife.h 3381 2010-11-15 12:33:00Z tforb $"
#define SEIFE_EXIT_FAILURE 1
#define SEIFE_CHECKERROR( EXPR , SUB, STR )\
if ( EXPR ) { fprintf(stderr, "ERROR (%s):\n %s\n", SUB, STR );\
exit(SEIFE_EXIT_FAILURE); }
/* seife filter types */
enum seife_filter_type {
IIRlpb, IIRhpb, IIRlp2, IIRhp2, IIRint, IIRhe1, IIRle1, IIRhe2, IIRle2,
IIRlp1, IIRhp1, IIRbp2
}; /* enum seife_filter_type */
/* global flags (to be initialized in cseife.c) */
extern struct seife_flags { int debug; } seife_global_flags;
/* structure to hold IIR filter coefficients */
struct seife_IIRcoef { double f0, f1, f2, g1, g2; };
/* enter debug mode */
void seife_debug_mode_on();
/* Butterworth lowpass (period t0, order o) */
void seife_lpb(double* x, int n, double dt, double t0, int o);
/* Butterworth highpass (period t0, order o) */
void seife_hpb(double* x, int n, double dt, double t0, int o);
/* 2nd order lowpass (period t0, damping h) */
void seife_lp2(double* x, int n, double dt, double t0, double h);
/* 2nd order highpass (period t0, damping h) */
void seife_hp2(double* x, int n, double dt, double t0, double h);
/* 2nd order bandpass (period t0, damping h) */
void seife_bp2(double* x, int n, double dt, double t0, double h);
/* 1st order lowpass (period t0) */
void seife_lp1(double* x, int n, double dt, double t0);
/* 1st order highpass (period t0) */
void seife_hp1(double* x, int n, double dt, double t0);
/* integration (time constant t0) */
void seife_int(double* x, int n, double dt, double t0);
/* 1st order highpass equalizer (former period t0s, new period t0) */
void seife_he1(double* x, int n, double dt, double t0s, double t0);
/* 1st order lowpass equalizer (former period t0s, new period t0) */
void seife_le1(double* x, int n, double dt, double t0s, double t0);
/* 2nd order highpass equalizer (former: period t0s and damping hs,
* new: period t0 and damping h)
*/
void seife_he2(double* x, int n, double dt,
double t0s, double hs, double t0, double h);
/* 2nd order lowpass equalizer (former: period t0s and damping hs,
* new: period t0 and damping h)
*/
void seife_le2(double* x, int n, double dt,
double t0s, double hs, double t0, double h);
/* detide with synthetic tides interpolated over ni samples */
void seife_tid(double* x, int n, double dt, int ni);
/* derivative (time constant t0) */
void seife_dif(double* x, int n, double dt, double t0);
/* set baseline to first value */
void seife_first(double* x, int n);
/* solve system of linear equations */
void seife_gauss(double *aik, int m, int n, double* rs, double* f);
/* apply IIR filter coefficients */
void seife_rekfl(double* x, double* y, int n, struct seife_IIRcoef c);
/* determine IIR filter coefficients */
struct seife_IIRcoef seife_rfk(enum seife_filter_type IIRtype,
double t0, double h, double t0s, double hs);
/* general Butterworth filter framework */
void seife_butterworth(double* x, int n, double dt,
enum seife_filter_type IIRtype,
double t0, int o);
#endif /* TF_CSEIFE_H_VERSION (includeguard) */
/* ----- END OF cseife.h ----- */
/*! \file cseife_deriv.c
* \brief calculate time derivative (implementation)
*
* ----------------------------------------------------------------------------
*
* $Id: cseife_deriv.c 3373 2010-11-14 13:29:20Z tforb $
* \date 28/06/2005
*
* calculate time derivative (implementation)
*
* Copyright 1984 by Erhard Wielandt
* This algorithm was part of seife.f. A current version of seife.f can be
* obtained from http://www.software-for-seismometry.de/
*
* REVISIONS and CHANGES
* - 28/06/2005 V1.0 Thomas Forbriger
*
* ============================================================================
*/
#define TF_CSEIFE_DERIV_C_VERSION \
"TF_CSEIFE_DERIV_C V1.0 "
#define TF_CSEIFE_DERIV_C_CVSID \
"$Id: cseife_deriv.c 3373 2010-11-14 13:29:20Z tforb $"
#include <cseife.h>
/* subs/seife_deriv.f -- translated by f2c (version 20000121).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
the code was derived through f2c, but modified thereafter
*/
/* derivative (time constant t0) */
void seife_dif(double* x, int n, double dt, double tau)
{
/* System generated locals */
int i__1;
/* Local variables */
static int j;
static double twodt;
/* derivative (non-recursive, non-causal, symmetric-difference) */
/* Parameter adjustments */
--x;
if (tau < 1.e-23) {
tau = 1.;
}
twodt = dt * 2. / tau;
i__1 = n - 2;
for (j = 1; j <= i__1; ++j) {
/* L1: */
x[j] = (x[j + 2] - x[j]) / twodt;
}
for (j = n - 1; j >= 2; --j) {
/* L2: */
x[j] = x[j - 1];
}
x[n] = x[n - 1];
} /* seife_dif */
/* ----- END OF cseife_deriv.c ----- */
/*! \file cseife_gauss.c
* \brief solve system of linear equations (implementation)
*
* ----------------------------------------------------------------------------
*
* $Id: cseife_gauss.c 3381 2010-11-15 12:33:00Z tforb $
* \date 28/06/2005
*
* solve system of linear equations (implementation)
*
* Copyright 1984 by Erhard Wielandt
* This algorithm was part of seife.f. A current version of seife.f can be
* obtained from http://www.software-for-seismometry.de/
*
* REVISIONS and CHANGES
* - 28/06/2005 V1.0 Thomas Forbriger
* - 15/11/2010 V1.1 do not use tfmacros.h
*
* ============================================================================
*/
#define TF_CSEIFE_GAUSS_C_VERSION \
"TF_CSEIFE_GAUSS_C V1.1"
#define TF_CSEIFE_GAUSS_C_CVSID \
"$Id: cseife_gauss.c 3381 2010-11-15 12:33:00Z tforb $"
#include <cseife.h>
#include <stdio.h>
#include <stdlib.h>
/* subs/seife_gauss.f -- translated by f2c (version 20000121).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
the code was derived through f2c, but modified thereafter
*/
void seife_gauss(double *aik, int m, int n, double* rs, double* f)
{
/* System generated locals */
int aik_dim1, aik_offset, i__1, i__2, i__3;
double d__1;
/* Local variables */
int imax[14];
double h__[15];
int j, k, l;
double q;
int index;
double aikmax;
SEIFE_CHECKERROR( m>13 , "seife_gauss", "matrix is too large" )
SEIFE_CHECKERROR( n>13 , "seife_gauss", "matrix is too large" )
/* solve linear equations */
/* Parameter adjustments */
--f;
--rs;
aik_dim1 = n;
aik_offset = 1 + aik_dim1 * 1;
aik -= aik_offset;
/* Function Body */
i__1 = m;
for (j = 1; j <= i__1; ++j) {
aikmax = 0.;
i__2 = m;
for (k = 1; k <= i__2; ++k) {
h__[k - 1] = aik[j + k * aik_dim1];
if ((d__1 = h__[k - 1], abs(d__1)) <= aikmax) {
goto L1402;
}
aikmax = (d__1 = h__[k - 1], abs(d__1));
index = k;
L1402:
;
}
h__[m] = rs[j];
i__2 = m;
for (k = 1; k <= i__2; ++k) {
q = aik[k + index * aik_dim1] / h__[index - 1];
i__3 = m;
for (l = 1; l <= i__3; ++l) {
/* L1404: */
aik[k + l * aik_dim1] -= q * h__[l - 1];
}
/* L1403: */
rs[k] -= q * h__[m];
}
i__2 = m;
for (k = 1; k <= i__2; ++k) {
/* L1405: */
aik[j + k * aik_dim1] = h__[k - 1];
}
rs[j] = h__[m];
/* L1401: */
imax[j - 1] = index;
}
i__1 = m;
for (j = 1; j <= i__1; ++j) {
index = imax[j - 1];
/* L1406: */
f[index] = rs[j] / aik[j + index * aik_dim1];
}
} /* seife_gauss */
/* ----- END OF cseife_gauss.c ----- */
/*! \file cseife_rekfl.c
* \brief apply IIR filter coefficients (implementation)
*
* ----------------------------------------------------------------------------
*
* $Id: cseife_rekfl.c 3373 2010-11-14 13:29:20Z tforb $
* \date 28/06/2005
*
* apply IIR filter coefficients (implementation)
*
* Copyright 1984 by Erhard Wielandt
* This algorithm was part of seife.f. A current version of seife.f can be
* obtained from http://www.software-for-seismometry.de/
*
* REVISIONS and CHANGES
* - 28/06/2005 V1.0 Thomas Forbriger
* - 11/07/2005 V1.1 support debug mode
*
* ============================================================================
*/
#define TF_CSEIFE_REKFL_C_VERSION \
"TF_CSEIFE_REKFL_C V1.1"
#define TF_CSEIFE_REKFL_C_CVSID \
"$Id: cseife_rekfl.c 3373 2010-11-14 13:29:20Z tforb $"
#include <stdio.h>
#include <cseife.h>
/* subs/seife_rekfl.f -- translated by f2c (version 20000121).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
the code was derived through f2c, but modified thereafter
*/
void seife_rekfl(double* x, double* y, int n, struct seife_IIRcoef coef)
{
/* System generated locals */
int i__1;
/* Local variables */
int j;
double xa, ya, xn, xaa, yaa;
/* report coefficients is requested */
if (seife_global_flags.debug == 1)
{
fprintf(stderr, "DEBUG (seife_rekfl):"
"%s=%-10.3g %s=%-10.3g %s=%-10.3g %s=%-10.3g %s=%-10.3g\n",
"f0",coef.f0, "f1",coef.f1, "f2",coef.f2, "g1",coef.g1,
"g2",coef.g2);
}
/* perform recursive filtering */
/* Parameter adjustments */
--y;
--x;
/* Function Body */
xa = 0.;
xaa = 0.;
ya = 0.;
yaa = 0.;
i__1 = n;
for (j = 1; j <= i__1; ++j) {
xn = x[j];
y[j] = coef.f0 * xn + coef.f1 * xa + coef.f2 * xaa + coef.g1 * ya +
coef.g2 * yaa;
xaa = xa;
xa = xn;
yaa = ya;
/* L1: */
ya = y[j];
}
} /* seife_rekfl */
/* ----- END OF cseife_rekfl.c ----- */
/*! \file cseife_rfk.c
* \brief calculate IIR filter coefficients (implementation)
*
* ----------------------------------------------------------------------------
*
* $Id: cseife_rfk.c 3373 2010-11-14 13:29:20Z tforb $
* \date 28/06/2005
*
* calculate IIR filter coefficients (implementation)
*
* Copyright 1984 by Erhard Wielandt
* This algorithm was part of seife.f. A current version of seife.f can be
* obtained from http://www.software-for-seismometry.de/
*
* REVISIONS and CHANGES
* - 28/06/2005 V1.0 Thomas Forbriger
*
* ============================================================================